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OBJECT 


•This  report  is  written  to  fulfill  the  requirements  of  Bureau  of  Naval 
Weapons  contract  SOv  65-0273-f  Item  4.  This  contract  was  issued  for 
development  of  a  digital  computer  program  ;o  predict  the  anti-icing 
requirements  for  engine  inlets  by  predicting  water  droplet  impinge¬ 
ment.  The  determination  of  the  design  requirements  is  made  by 
establishing  the  potential  flow  field  around  the  subject  surface  and 
then  calculating  the  trajectory  of  the  water  droplet  as  it  approaches 
and  impinges  on  the  surface. 


During  the  execution  of  this  contract,  any  new  discoveries  or  inventions 
(subject  inventions)  were  to  be  reported  to  the  Government. 
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SUMMARY 


2.0 


2.1  Development  Summary 

With  the  signing  of  the  development  contract  on  March  9)  1965,  the 
Contractor  started  work  on  the  two  "basic  digital  computer  programs  to 
provide  an  accurate  method  of  water  droplet  impingement.  The  potential 
flow  computer  program  was  updated  and  refined  "by  the  addition  of  vari¬ 
able  distance  incrementation,  improved  plotting  and  table  preparation 
routines  and  more  output. 

The  water  droplet  trajectory  computer  program  was  adapted  to  the 
improved  potential  flow  program  and  expanded  to  include  methods  to  cal¬ 
culate  trajectories  that  included  smaller  droplets,  droplet  reversal, 
gravity  and  buoyancy  effects  and  random  input.  Provisions  for  summing 
the  results  from  a  given  droplet  size  distribution  and  for  determining 
the  tangent  droplet  trajectories  were  also  included.  Both  computer 
programs  were  verified  by  comparison  with  theoretical  and  NACA  analytical 
and  test  results. 

In  compliance  with  contract  Section  V  (FATENT  RIGHTS),  the  Contractor 
hereby  certifies  that  there  are  no  Subject  Inventions  relative  to  the 
fulfillment  of  this  contract. 

2.2  Descriptive  Summary 

The  Boeing  Company  under  contract  to  the  Bureau  of  Naval  Weapons  has 
developed  a  digital  computer  program  to  predict  the  anti-icing  require¬ 
ments  of  engine  inlets  by  calculating  water  droplet  trajectories  toward 
the  particular  inlet  through  a  potential  flow  field.  The  Incompressible 
potential  flow  field  is  determined  by  a  separate  computer  program  using 
known  boundary  conditions.  The  potential  flow  field  results  are  combined 
with  the  droplet  and  flow  physical  data  for  input  to  a  droplet  trajectory 
computer  program.  This  water  droplet  trajectory  program  is  based  on 
balancing  the  dreg,  buoyancy  and  weight  forces  on  the  water  droplet  with 
the  rate  of  change  of  the  water  droplet  momentum. 

The  water  droplet  computer  program  is  extremely  versatile  because  the 
actual  design  shape  is  used  in  the  study  and  therefore  the  results  do  not 
need  to  be  extrapolated. 

Verification  of  the  computer  programs  was  accomplished  by  comparisons 
with  theoretical  and  NACA  test  and  analytical  data. 
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3.0  CONCLUSIONS 

Based  on  historical  and  current  development  experience  as  shown  by 

the  references  and  this  report,  it  is  concluded  that:  • 

1.  This  method  of  determining  water  impingement  on  cylindrical  and  plane 
objects  is  accurate. efficient  and  rapid. 

2.  The  water  droplet  trajectory  computer  program  can  be  used  to  predict 
the  trajectories  of  large  and  small  water  droplets  starting  at 
arbitrary  conditions. 

3.  The  potential  flow  computer  program  is  a  useful  tool  to  describe 
complex  flow  fields  to  assist  in  surface  and  duct  pressure  distri¬ 
bution  studies. 

k.  The  developed  computer  program  does  not  contain  any  new  inventions 
or  discoveries  that  could  be  classed  as  a  Subject  Invention  (contract 
section  V) . 
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4,0  RECOMMENDATIONS 

Based  on  the  results  and  conclusion  described  in  this  document, 

it  is  recommended  that? 

1.  This  method  of  water  impingement  prediction  be  adopted  as  a 
standard  for  anti-icing  requirements  determination, 

2.  Development  of  the  water  droplet  trajectory  program  be  continued 
to  expand  its  utility,  simplify  its  usage  and  increase  its 
accuracy  especially  for  heavier,  larger  particles, 

3.  The  potential  flow  computer  program  be  expanded  to  Include 
compressibility  effects. 
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5.0  DESCRIPTION 


A  preliminary  water  droplet  trajectory  calculation  computer  program  was 
refined,  expanded  and  verified  during  the  completion  of  Item  1  of  the 
subject  contract.  Specific  effort  was  devoted  to  accomplishing  the 
following, 

(a)  The  potential  flow  field  computer  program  was  updated  and  refined 
for  increased  utility  and  versatility.  Variable  distance  incre¬ 
mentation  in  both  the  horizontal  and  vertical  (or  radial)  directions 
was  provided  to  more  adequately  define  a  potential  flow  field  by 
using  fine  increments  in  areas  of  special  interest.  The  potential 
flow  program  output  was  also  increased  to  include  velocity,  velocity 
ratio  and  pressure  coefficients  at  coordinate  points  on  selected 
streamlines . 

(b)  The  water  droplet  trajectory  computer  program's  capabilities  were 
expanded  to  permit  the  calculation  of  smaller  droplet  trajectories 
and  droplet  direction  reversal.  These  changes  were  accomplished 

by  making  the  water  droplet  horizontal  distance  increments  independent 
of  those  selected  for  the  potential  flow  grid  and  by  using  an  addi¬ 
tional  set  of  equations  for  vertical  travel.  For  droplet  direction 
change  calculations,  the  droplet  primary  direction  of  travel  is 
changed  whenever  the  perpendicular  distance  traveled  is  greater 
than  an  element  in  that  direction.  The  droplet  direction  is 
reversed  whenever  the  droplet  inertia  is  insufficient  to  carry  it 
across  an  increment  and  the  perpendicular  velocity  component  is 
small. 

(c)  The  water  droplet  trajectory  computer  program  was  modified  to  pro¬ 
vide  for  variation  in  the  water  droplet's  initial  state  with  respect 
to  the  free  stream  flow  (random  input)  and  for  the  effects  of 
gravity  and  buoyancy  on  the  droplet  trajectories.  The  random  input 
feature  permits  starting  the  droplet  trajectory  from  any  point  in 
the  potential  flow  field  with  any  velocity  components.  The  inclu¬ 
sion  of  gravity  and  buoyancy  effects  destroys  the  symmetry  of  the 
cylindrical  flow  field  trajectories  and  requires  calculations  with 
negative  ordinate  values.  The  capability  of  calculating  below  the 
centerline  has  been  included. 

(d)  The  water  droplet  trajectory  computer  program  was  verified  by 
comparison  of  the  program's  impingement  predictions  with  NACA 
analytical  and  test  results.  Analytical  comparisons  were  made 
against  the  NACA  cylinder,  sphere  and  airfoil  data  and  comparisons 
with  the  test  data  were  made  for  the  sphere,  airfoil  and  supersonic 
engine  inlet. 
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6.0  MFffBDD 

This  analytical  method  of  water  droplet  impingement  prediction  utilizes 
a  digital  computer  and  a  finite  difference  solution  of  the  trajectory- 
equations  of  particles  in  an  incompressible  potential  flow  field.  The 
potential  flow  field  is  determined  from  known  boundary  stream  function 
values  by  using  the  relaxation  technique  for  solving  the  stream  function 
equations  derived  from  continuity  and  irrotationality. 

The  results  of  the  potential  flow  computer  program  are  used  with  addi¬ 
tional  water  droplet  and  flow  physical  data  as  input  to  the  water  droplet 
trajectory  computer  program.  The  water  droplet  trajectory  program 
calculations  result  in  an  accurate  determination  of  water  impingement 
intensity  on  a  desired  surface. 

Flow  charts  and  listing  of  the  computer  programs  are  presented  in 
appendices  1  through  4. 

6,1.1  Potential  Flow  Field  Analysis 

Since  the  computer  program  solves  the  stream  function  equations  to  obtain 
the  potential  flow  field  from  boundary  stream  function  values,  these 
boundary  values  must  be  determined  analytically.  The  generalized  flow 
field  for  an  engine  inlet  is  shown  on  page  12.  Using  this  figure  with 
the  y- values  as  radial  lenghts,  we  can  develop  the  boundary  values  as 
follows • 

For  two  dimensional  (plane)  incompressible  steady  fluid  flow,  the  velocity  . 
component  in  any  direction  is  found  by  differentiating  the  stream  function 
at  right  angles  to  that  direction.  Assuming  the  positive  flow  direction 
as  upward  (increasing  y-values)  and  to  the  right  (increasing  x-values) 
the  velocity  components  are 


/V^=  _ 

(1) 

AT  =  - 

jy' 

(2) 

J  x 

In  the  case  of  axi ally-symmetric  flow,  such  as  that  of  an  engine  inlet 


the  velocity  components  are 

SV^i  r.  1  fll  = 

(3) 

y 

o>  y 

/Vy  =  /^"radial  =  - 

6Y 

(*) 

y  <9X 
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The  airflow  toward  the  engine  inlet  is  assumed  to  be  from  left  to  right 
with  a  value  of  U  &  at  the  left  boundary.  At  this  boundary  the  flow 
velocity  is  uniform  and  horizontal.  Therefore, 


• 

/Vyi  =  /Vy  left  =  0 

(5) 

'Vxl  -  left  -  c,'^=UL  =  Uo5=kui? 

y  6  y 

(6) 

then 

J*  ydy 

(7) 

and 

y  =  h  u  ‘co  y  2/2  +  c 

(8) 

or 

y/l  =  k  u'co  y i  %  +  C 

(9) 

To  evaluate  the  constants  (  k  and  c)  the  following  boundary  conditions 
are  used. 

(a) 

'V/i:=0aty].  =  0 

(b) 

Yl  =  1  at  yj_  =  Rco  (stream  tube  radius) 

Using 

(a)  c^  is  evaluated  as  zero. 

2 

Using  (b)  k  is  evaluated  as  2/  (  Uo>  R  co  )  and  equation  (9)  can  be 
rewritten  as 

Ti  =  y  2/  2 

yl  ‘  RCO 

(10) 

At  the  radial  extremity  of  the  control  volume  (y  =  y  max) 

u 

S 

ii 

1 

Since  ymaoc  a  constant  and  assumed  to  be  the  ordinate  of  a  horizontal 
stream  line,  y/  max  a  constant  for  all  values  of  x. 

The  arbitrary  selection  of  the  stream  function  value  of  one  (l.O)  for  the 
value  at  the  stream  tube  radius  on  the  left  boundary  establishes  this 
value  everywhere  on  the  stream  tube  boundary  including  the  cowl  surface. 
Similarly,  selection  of  the  zero  stream  function  value  on  the  centerline 
establishes  this  value  on  the  centerline  and  on  the  centerbody  surface. 
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^  At  the  right  boundary,  uniform  horizontal  flow  is  assumed  to  exist 

within  the  inlet  and  also  outside  the  cowl.  These  flows  can  (and 
generally  do)  have  different  velocities.  With  these  assumptions 
the  following  basic  equation  can  be  developed  and  evaluated. 

7  =  k  U  y2/2  +  c.  '  '  (II) 

\ 

Inside  the  engine  inlet 

7=  0  at  yr  =  y  ocb  so  c2  =  -  kg  u  y2ocb/2- 

2 

y  =  1.0  at  yr  =  yic  so  k2  =  2/(u  (yic  -  yj^)). 

Then 

/  2  2  ,  2  2  .  .  . 

Y ri=  (  yr  ■  y0cb)  /  (  y  ic  “  y  ocb)  (12) 


This  equation  for  the  stream  function  is  valid  between  the  engine  inlet 
centerbody  and  cowl  at  the  right  side  of  the  control  volume. 

Outside  the  engine  inlet  equation  (ll)  is  evaluated  with  these  condi¬ 
tions  . 


y- 

-•  1.0  at  yr 

=  yoc 

r- 

.  iy  _ 
max-  p 

at  yr 

=  y 

max 

R-5 

2  , 

Then 

c3  = 

:  1  -  1 

^  ^oc  / 

2 

k3  = 

=  (2/  U  )  ( 

y2  /  r2 

y  max'  53 

-1)/ 

^^max 

2 

-  yoc) 

Y 

ro  ~ 

i  +  (yr2  - 

y0c) 

(  2  /  2 
'■Ymax  co 

-1)/ 

(y  2~ 
wmax 

y0c) 

(13) 


This  equation  is  used  to  evaluate  the  stream  function  values  on  the  right 
control  volume  boundary  in  the  constant  velocity  region  outside  the 
engine  inlet. 

For  a  particular  engine  inlet  configuration  and  engine  power  setting 
(airflow  requirement),  the  input  left  and  right  boundary  stream  function 
values  for  the  potential  flow  program  are  determined  from  equations  (10), 
(12)  and  (13).  The  lower,  upper  and  model  surface  stream  function 
boundary  values  are  arbitrarily  assigned  at  zero,  maximum  and  one, 
respectively. 
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6.1.2  Potential  Flow  Computer  Program 

Incompressible,  inviscous  and  irrotational  fluid  flow  problems  in  plane 
and  axi ally-symmetric  flow  fields  can  be  solved  with  the  assistance  of 
the  IBM  709V  FORTRAN  digital  computer  program  described  herein. 

Input  data  for  the  program  includes  a  geometric  description  of  the  flow 
field  boundaries  and  stream  function  values  for  the  boundaries.  Program 
results  include  a  complete  description  of  the  potential  flow  field  by 
enumerating  stream  function  values  at  mesh  intersections  and  streamline 
coordinates,  as  well  as  fluid  velocities  and  pressure  coefficients  for 
selected  streamlines. 

The  stream  function,  psi  (y),  for  incompressible,  steady,  fluid  flow  is 
defined  as  the  volume  of  fluid  passing  between  a  base  streamline  and  the 
point  in  question.  For  a  plane  flow  field  the  field  depth  ie  considered 
to  be  one  unit  and  for  an  axially  symmetric  flow  field,  the  depth  varies 
with  the  radius. 

The  velocity  component  in  any  direction  is  the  differential  of  the  stream 
function  at  right  angles  to  that  direction.  The  usual  convention  of 
positive  flow  upward  and  to  the  right  results  in  the  following  derivatives 
for  velocity  components. 

Field  Plane  Axially  Symmetric 

Horizontal  velocity,  Vx  o^/cJy  (l/y)  (  f)  Yf  £  y) 

Vertical  velocity,  Vy  -  dYf  x  (-l/y)  (  d  ^7  a x) 


Limiting  the  flow  to  be  irrotational  introduces  the  following  equation: 

aV*/ay  -o 

^y/$X2+  ^y/dt  »  0  (14) 

for  plane  flow.  For  axially  symmetric  flow  this  limitation  yields 
■>/ax  (-Vy  >r/d>.  )  -  VaY(W  )  =  0 
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Equations  (l4)  and  (15)  define  the  stream  function  for  incompressible, 
steady,  inviscous  and  irrotational  fluid  flow  in  plane  and  axially 
symmetric  flow  fields.  Numerical  evaluation  of  these  equations  by  the 
relaxation  technique  is  the  primary  function  of  the  potential  flov 
computer  program. 


The  potential  flow  computer  program  is  divided  into  three  distinct  sections: 
potential  flow  field  solution;  derivative  calculation  and  contour  mapping. 
The  potential  flow  field  solution  is  accomplished  by  a  finite  difference 
approximation  of  equations  (l4)  and  (15)»  This  solution  is  based  on  the 
use  of  a  rectangular  grid  with  the  origin  of  coordinates  temporarily  at 
the  point  being  considered.  For  a  typical  mesh  point,  fo  is  the  stream 


An  example  of  a  mesh  point  coordinate  system  showing  variable  mesh  spacing 
and  coordinate  orientation  is  illustrated  by  this  sketch.  The  average 
first  partial  derivaties  of  "VC  with  respect  to  x  from  the  origin  to 
points  1  and  3  are: 

dr/)x)0>1  'Ar/ix>o,i ‘  -Xo)A  .  (16) 

a  r/JX>3,0  =  Ar/^3,0  "  Wo  -r3)/  0  (IT) 


The  average  second  partial  derivative  of  vith  respect  to  x  to  be  used 
between  the  distance  a/2  to  the  right  of  the  origin  to  c/2  to  the  left 
of  the  origin  is 


JBT£FA FS/V£?\ 
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<>2  r/sy?)0  =  (ay/d  x)  0j  x  -  dy/dx)  3j  0)/ (a  ♦  e)  /  2 

=  2(V'i  -y/Q)/(a  (a+c))  -  2  (TQ  ->3)/  (  c  (afc))  (l8) 

By  similar  analysis  in  the  vertical  direction 
2 

6  y/  d y  \  =  2  (y  2-y0 )/  o>o>+4))  -  s(y0  -r4)/(d(b  ■fd»  («) 

The  first  partial  derivative  of  y;  with, respect  to  y  at  the  origin  is 
approximated  by  a  linear  interpolation  of  the  average  first  partial 
derivaties  in  the  adjacent  intervals  above  and  below  the  origin. 


d  Y/d  y)8.  -  mr/d  y)  0,  2  +  »  ( ar/d  y)  ki0V  <b  +  d> 

-  a  (r2-y0)  /  0>0*a))  + *  {fo  -Yk>  /  W^4))  (2°) ' 


Substitution  of  equations  (18),  (19)  and  (20)  into  equations  (14)  and  (15) 
results  in  the  following  finite  difference  approximations  for  plane  and 
axially  symmetric  flow,  respectively.. 


gyi  +  g^a  + SY 3  +  ZY\ 

a(a+c)  b  (b+d)  c  (a  +c)  d(b  +d) 

2  2  2  2 
a"(a  +c)  +  b  (b+d)  x  c(a  +c)  +  d  (b  +  d) 


(21) 


To  = 


Yx 


(2-d/y)y/  +  2y"3+  (2  +  b/y)’A\ 


a(a+c) 


b  (b+d) 
(2-d/y) 


c  (a+c)  d(  b  +  d) 

2  (2  +b/y) 


(22) 


a(a+c) 


b(b+d) 


c(a+c) 


d(b  +d) 


In  the  computer  program  equations  (22)  is  reduced  to  equation  (21)  for 
plane  flow  by  making  d/y  and  b/y  equal  zero. 


The  computer  program  constructs  a  finite  difference  solution  at  each  mesh 
point  and  sweeps  through  them  in  a  definite  order.  Repetitive  calculations 
(sweeps)  result  in  the  stream  function  values  converging  to  values  consistent 
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with  the  specified  boundary  values*  To  improve  convergence  over-relaxation 
is  used.  Shat  is,  the  calculated  change  in  the  stream  function  value  at 
each  point  is  multiplied  by  a  relaxation  factor  which  is  greater  than  one. 
This  result  is  added  to  the  previous  value  and  the  sum  is  used  for  the  nev 
stream  function  value* 


Sweeping  is  done  in  alternate  directions,  first  by  columns  from  lower  left 
to  upper  right,  and  then,  also  by  columns,  from  upper  right  to  lower  left. 
(This  improves  the  rate  of  convergence  somewhat  for  cases  with  normal 
derivative  enforcement  on  the  boundaries) .  A  series  of  to  sweeps  (20  pairs 
of  alternate  direction  sweeps)  is  carried  out  to  establish  the  pattern  of 
)  convergence;  then  the  point  with  maximum  error  is  located  and  followed 

through  ten  more  sweeps.  The  change  in  the  indicated  error  over  this  last 
sequence  of  eleven  sweeps  is  used  to  define  the  ratio  for  a  geometric  series. 
The  ratio  is  used  to  extrapolate  the  migration  of  the  field  over  the  sequence 
of  eleven  sweeps  to  an  indicated  limit,  subject  to  a  hedge  factor  and  a 
maximum  limit  on  the  indicated  extrapolation.  This  procedure  has  been  found 
to  accelerate  the  convergence  of  well  behaved  configurations  by  about  30$ 
and  the  convergence  of  poorly  behaved  configurations  by  a  factor  of  three. 

The  second  pert  of  the  computer  program  operates  on  the  potential  flow 
field. obtained  above  and  calculates  the  partial  derivatives  of  with 
respect  to  x  and  y  at  all  mes'n  points.  These  derivatives  are  obtained  by 
linear  interpolation  using  the  values  at  the  desired  point  and  at  the 
adjacent  points  on  each  side  at  the  same  ordinate  or  abscissa  value. 
Derivatives  at  points  on  boundaries  are  extrapolated  linearly  from  the  two 
points  closest  to  the  boundary  at  the  same  ordinate  or  abscissa-yalup . 

The  absolute  velocity  at  each  mesh  point  is  calculated  as  the  square  root 
of  the  sum  of  the  squares  of  the  x  and  y  partial  derivatives  of  at  the 
point.  For  axially  symmetric  flow  fields  this  calculation  includes  division 
by  the  y  value  (radius)  tc  obtain  the  true  velocity. 

Contour  mapping  of  selected  stream  lines  or  constant  partial  derivatives 
is  accomplished  in  the  third  program  section.  The  ordinate  values  of  stream 
lines  (stream  function  isolines)  or  contours  of  constant  partial  derivative 
values  are  determined  at  t.is  input  abscissa  values  from  the  potential  flow 
field  and/or  either  of  the  derivative  fields  by  linear  interpolation.  The 
flow  velocity  at  the  stream  line  coordinates  is  also  determined  by  inter- 
>  polation. 

The  pressure  coefficient  is  calculated  at  each  point  on  the  contour  lines. 

The  pressure  coefficient  is  calculated  as  1  -(local  velocity/ ffee  stream 
velocity) 2. 
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6.1.3  Potential  Flow  Computer  Program  Usage 


Input  for  the  potential  flov  computer  program  is  prepared  on  six  field, 
72  column  floating  decimal  data  cards  as  illustrated  on  the  following 
Potential  Plow  Program  Input  Format.  The  individual  input  quantities 
are  described  as  follows. 


TITLE  Identification  of  the  problem  -  used  as  a  heading 

in  the  printout. 

AXYM  Flow  field  type  key  -  zero  for  a  plane  flow  field  and 

one  for  the  axially  symmetric  flow  field. 

CVGS  Relaxation  convergence  tolerance  -  when  the  maximum 

difference  between  two  successive  sweep!  s  v  PSI  value  .is 
less  than  this  number,  the  flow  field  iteration  will 
stop.  This  value  should  be  greater  than  3  x  10 
due  to  machine  limitations . 

NSWPS  Number  of  sweeps  -  maximum  number  of  sweeps  permitted. 

The  iteration  will  stop  after  this  number  of  sweeps. 

PUNCH  Potential  flow  field  table  punch  key  -  table  will  be 

punched  if  this  value  is  one. 

NDERIV  Normal  derivative  key  -  one  to  force  the  normal  deriva¬ 

tives  at  boundaries  to  zero.  Normally  input  as  zero 
and  yo  and  II  are  ignored. 

YO  Y-value  of  the  lower  boundary  on  which  the  normal  deriva¬ 

tives  are  forced  to  zero. 

71  Y-value  of  the  upper  boundary  on  which  the  normal  deriva¬ 

tives  are  forced  to  zero. 

*XSCL  Scale  factor  in  the  X  direction  for  plot  of  stream  lines 

equal  to  the  number  of  inches  which  correspond  to  a  unit 
change  in  X. 


*YSCL 

Scale 

*XL1 

Lover 

*XL2 

Upper 

*711 

Lower 

*YL2 

Upper 

*  XSCLA,  YSCLA,  XL1A 
for  a  second  plot, 
be  blank. 

REVLTR: 


factor  in  the  Y  direction. 

limit  in  the  X  direction  of  area  to  be  plotted, 

limit  in  the  X  direction  of  area  to  be  plotted. 

limit  in  the  Y  direction  of  area  to  be  plotted. 

limit  in  the  Y  direction  of  the  area  to  be  plotted. 

XL2A,  YL1A  and  YL2A  are  the  corresponding  parameters 
If  only  one  plot  is  desired  the  sixth  card  should 
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Yl,  Y2,  Y3, 

...  Ym 

First,  second,  third  ...  m'th  selected  ordinate 
mesh  line  values. 

Xa,  Xb,  Xc, 

***  \/2 

First,  second,  third  ...  n/2'th  abscissa  boundary 
values  at  a  given  ordinate  value. 

7a,  Tb,  To, 

’••Yn/2 

Stream  function  boundary  values  at  Xa,  Xb,  Xc,  ... 
Xjj^g  or  at  Ya,  Yb,  Yc,  ...  Y x\.j2* 

:  n 

The  number  of  boundary  coordinates  and  stream 
function  values  along  that  particular  mesh  line. 

XI,  X2,  X3, 

•  •  •  Xm 

First,  second,  third  ...  m*th  selected  abscissa 
mesh  line  values. 

Ya,  Yb,  Yc, 

•••  Yn/2 

First,  second,  third  ...  n/2‘th  ordinate  boundary 
value  at  a  given  abscissa  value. 

N3L 

Number  of  stream  lines  of  a  type  (see  CSCH). 

CSCH 

Determines  whether  stream  lines  connect  points  of 
equal  stream  function  value  (CSCH  =2),  derivative 
of  PSI  with  respect  to  X  (CSCH  *  3)  or  derivative  of 
PSI  with  respect  to  Y  (CSCH  =  4). 

NNLNS 

Number  of  types  of  SLINES.  NSL,  CSCH,  and  the  SLINES 
must  be  input  NNLNS  times. 

SLINES 

Values  of  PSI  (  or  derivatives  of  PSI  with  respect 
to  X  or  Y)  through  which  lines  will  be  drawn. 
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FLOATING  DECIMAL  DATA  FOR  7094 


Output  from  the  potential  flow  computer  program  is  self  explanatory 
and  is  printed  in  the  following  order. 

Input  Input  is  printed  out  to  provide  a  permanent 

record. 

Sweep  history  Record  of  convergence  is  printed  out  with 

residue  at  point  of  largest  deviation. 

Potential  flow  field  solution  -  The  stream  function  values,  deriva¬ 
tives  and  flow  velocity  are  printed  out  at 
each  mesh  point. 

Contour  lines  Stream  line  or  constant  partial  derivative 

coordinate  values,  velocities  and  pressure 
coefficients  at  these  points  are  printed. 

On  request,  the  program  will  also  provide  a  plotter  tape  for  a  Gerber 
plotter  of  the  contour  lines.  The  potential  flow  field  is  punched 
on  data  cards  in  table  form  for  use  in  the  water  droplet  trajectory 
computer  program. 


The  equations  of  motion  for  a  water  droplet  in  air  are  obtained  by 
a  summr.tion  of  forces  on  the  drop.  Ignoring  the  gravitational  force, 
the  drag  force  on  the  water  drop  must  be  equal  and  opposite  to  the 
rate  of  change  of  the  water  drop  momentum.  The  drag  force  in  the 
x-direction  is 

Dx  =  <13  (  P*.  /Vfel  /S)  (  TT<?)  (ATX  -  1  (23) 

where  rel  =  J\T  ~/\Ty  . 

The  rate  of  change  of  drop  momentum  in  the  x-direction  is 

Fx  =  mass  x  acceleration  =  (4  Tf  /3)  dyi/^/dt.  (24) 
The  drop  Reynolds  number  is 

R  =  (2  a  A  ^ vel^U(  (25) 

then  Dx  =  CD  R  (/}/£  *  yi fa)/  ^  '  (26) 

and  Fx  =  Dx  (27) 

(4  Tf  a3  A/3)  d/VVdt  =  (CD  rA)  Tf  a A  W'x  -/y£v)  (28) 

or  (2  a2 u  J9A-  )  d/lTxw/dt  =  (  U  CjjR/24)  (/V^  "  /Txw)  (29) 
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where  /V^o,  =  0/4rxv>  Srx  - 

i'rxv='*W^l 


d/tW  dt  =  (UCdR/24k)  (/it  -/y^) 

(30) 

where  K  =  2  ^  f  jl9jA 

(31) 

and  since  dt  =  dx/zV^U 

(32) 

'V'xw  d<  =  (CdR/2^)  (/T^  -/Txv)  to. 

(33) 

Similarily  in  the  y-direction 

d/^/dt  =  UCjjRMk  (/i/j.  -/y'yW) 

(34) 

Equation  (33)  can  be  integrated  using  mean  values  of  velocity  over 
an  interval. 

(/^xw)n+l 

Xn+1 

J/VZx  d/l^w  =  (CjjRMk)^^  A \ 

/dx 

(35) 

(/Vxw)n 

Xn 

JWlcw)  n+1- 

/2 

=  (CDK/Ste)^^  (/rx  -/fe,)mean(xn+i  -  xn)  (36) 


Equation  (3*0  can  be  integrated  in  a  similar  manner. 


REVLTR: 


(/^yv^  n+1 


'^yv)n 


(UCoR/SliKj^^ 


(-n  -/n*> 


yw/mean 


/a 


Jn+1 

dt 


(37) 


^n 

4BF£J&'fA/£? 
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C'Vwl  -  {ArJ  n  ■  ^  -^)JV a-tj  (38) 

From  equation  (32)  is  obtained 


and  similarily 

4+1  "  4  +■ «  |^v)n+l  +  (^n  (4+1  ‘  4)/  2-  W 

By  similar  analysis  using  mean  values  of  velocity  over  an  interval  in 

the  y-direction  the  following  equations  are  derived. 

> 

[^yv^n+l  “  (/^yv).  ]/  2  )mean  “/^v^mean  (yn+l  "  Yn)  ^1) 

4 

(/^xw)n+l  "  (/l/xv)n  =  (UCpR/24K)mean(/2^c  -/KxwJjueanCtn+i  "  tn)  (42) 

^n+1  ~  2(^n+l  "  ^n)/u( (/^/yv)n+l  +  (/^/yv)n)  (^3) 

Xn+1  =  Xn  +  U  |^xw)n+l  +  (/^xw)n]  ^n+1  **  ^n)/2  (^0 

Equations  (36),  (38),  (39)  and  (40)  are  the  basic  equations  of  the  water 
droplet  trajectory  computer  program  for  droplet  travel  that  is  predominately 
in  the  horizontal  direction.  For  droplet  travel  in  the  vertical  direction 
equations  (4l)  through  (44)  are  used.  The  use  of  equations  for  travel  either 
horizontally  or  vertically  permits  the  calculation  of  droplet  trajectories 
in  areas  where  a  fixed  interval  is  not  traversed  in  the  selected  direction. 

The  inclusion  of  gravitational  and  buoyancy  effects  on  the  water  droplet 
trajectories  results  in  the  following  modifications  to  the  trajectory 
equations.  Equation  (38)  contains  the  additional  term  -(g/y^yU) 

(tn+i  -  tn)  on  the  right  side.  Equation  (4l)  must  include  the  term 

"(g/Z^wU2)  K  "fa)  Un+1  ~  yn)  on  the  riSht  side« 
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6.2.2  Water  Droplet  Trajectory  Computer  Program 

The  water  droplet  trajectory  computer  program  traces  the  trajectories 
of  individual  drops  across  finite  incremental  distances  through  the 
incompressible  potential  flow  field  until  they  impinge  on  a  surface 
or  pass  out  of  the  area  of  interest.  This  program  includes  ail” 
velocity  components  and  trajectory  routines  with  appropriate  tests 
as  shewn  on  the  following  abreviated  flow  chart  and  on  the  detail 
flow  chart  and  listing  of  the  appendix. 


A  droplet  size  distribution  impingement  intensity  summation  routine 
was  added  to  the  water  droplet  trajectory  program  as  an  option. 

This  routine  was  added  to  provide  a  more  accurate  comparison  of  the 
program  re  stilts  with  experimental  data  in  artificial  and  natural 
icing  conditions.  The  Langmuir  "D"  droplet  size  distribution  of 
reference  1  was  selected  and  incorporated  in  the  program.  Other 
droplet  distributions  could  be  added  with  small  program  chang, es. 

The  Langmuir  "D"  distribution  is  as  follows. 


Liquid  water  in  size  group, 
percent _ 


Average  group  drop  radius/volume 
mean  drop  radius _ 


5 

0.31 

10 

0.52 

20 

0.71 

— f 

30 

1.00 

20 

1.37 

10 

1.74 

5 

2.22 
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6.2.3  Water  Droplet  Trajectory  Computer  Program  Usage 


Input  for  the  water  droplet  trajectory  computer  program  consists  of 
geometric  and  potential  flow  descriptive  tables  and  data  cards.  The 
tables  are  read  first  in  numerical  order  from  one  through  nine. 

Table  arrangement  is  shown  on  the  following  pahle  Input  Format . 
Symbols  used  in  the  Table  Input  Format  are  described  as  follows. 


YIvC 


MG 


Number  of  ordinate  values  in  the  table  plus 
one,  (KGTB:  Y2i0  is  2.  for  a  single-  function 
table . ) 

Number  of  abscissa  values  in  the  table  plus 
one* 


Table  No. 


Ig  •••  Yn 


•  *  « 


Index  number  of  table. 

First,  second  ...  n'th  ordinate  values 
First,  second  ...  n*th  abscissa  values 


*(*i>  *j> 


Dependent  variable  value  at  (Xi,  Y^)  where  i  has 
values  1,  2,  ...  m  and  j  has  values  1,  2,  ...  n. 


The  physical  orientation  of  the  individual  tables  is  shown  on  the  Water 
Droplet  Trajectory  Computer  Program  Table  Key.  A  written  description 
of  these  tables  is  as  follows. 


▼ 
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Water  Droplet  Trajectory  Computer 
Program  Table  Key 


Engine 


Table  4 
Table  5 


—  4* 


YBl  Table  8* 
YB2  Table  6** 
YB3  Table  7** 


*  Table  9  is  generated  from  Table  8  from  the  centerbody  highlight  aft. 

**  Tables  1  and  2  are  generated  from  Tables  6  and  7«  respectively, 
from  the  cowl  highlight  aft* 


DATE  i 


8  CONTRACT  HO. 
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Table  Nmsber  Variable  Description 

1  Lower  cowl  surface  distance  from  highlight  versus 

horizontal  distance  to  right  of  highlight,  A 
dummy  table  is  input  and  the  actual  values  are 
generated  internally  from  Table  6  values. 

(SUP  vs.  X-XEL) 

2  Upper  cowl  surface  distance  from  highlight  versus 

horizontal  distance  to  right  of  highlight.  A 
dummy  table  is  input  and  the  actual  values  are 
generated  internally  from  Table  7  values. 

(SUP  vs.  X-XHL). 

3  Stream  function  values  versus  X  and  Y.  This  is 

the  only  double  function  table  used  and  it  is 
generated  by  ;he  potential  flow  program.  Hie  X  and 
Y  values  used  are  the  mesh  points  (  y*  vs.  X,  Y) 

k  First  boundary  value  of  X  after  left  hand  boundary 

at,  given  Y  value  (XB2  vs.  Y). 

5  Second  boundary  value  of  X  after  left  hand  boundary 

at  given  Y  value.  (XB3  vs.  y) 

6  First  boundary  value  of  Y  above  lower  boundary  at 

given  X  value.  (YB2  vs.  X). 

7  Second  boundary  value  of  Y  above  lower  boundary 

at  given  X  value  (YB3  vs.  X) 

8  lowest  boundary  value  of  Y  at  given  X  value.  This 

is  normally  an  extension  of  the  engine  centerline 
and  the  centerbody  surface  (YB1  vs.  X) 

9  Centerbody  surface  distance  from  highlight  versus 

horizontal  distance  to  right  of  highlight. 

A  dummy  table  is  input  and  the  actual  values  are 
generated  from  table  8  values.  (SUP  vs.  X-XKL). 

After  the  tables  are  read  into  the  computer,  the  input  data  is  called. 

These  data  are  arranged  on  six  field  72-column  floating  decimal  cards 
as  shown  on  the  following  Water  Droplet  Trajectory  Computer  Program 
Input  Format.  The  Individual  input  quantities  are  described  as  follows. 
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Tape  No. 

X  scale  1 

Y  scale  1 
X  set  1 

Y  set  1 
X  max  2 

Y  max  2 

X  scale  2 

Y  scale  2 
X  set  2 

Y  set  2 
LWC 

a 

Uoo 

*4 

fa 
/<a 
Xmsx 
Yraax 
Ydrop 
A  Y3rop 

XL 

XKL1 


Problem  title  -  used  as  a  heading  on  printout 
and  title  on  machine  jlot. 

Index  number  of  plotting  tape  unit  (usually  6) 
X-scale  factor  for  first  curve 
Y-scale  factor  for  first  curve 
Minimum  X-value  of  first  curve-inches 
Minimum  Y-value  of  first  curve-inches 
Maximum  X-value  of  second  curve-inches 
Maximum  Y-value  of  second  curve  -  inches 
X-scale  factor  for  second  curve 
Y-scale  factor  for  second  curve 
Minimum  X-value  of  second  cuive -inches 
Minimum  Y-value  of  second  curve-inches 
Liquid  water  content  of  air-lb/cu.  inch 
Droplet  radius  -  ft. 

Actual  remote  free  stream  velocity  -  ft/sec. 

Potential  flow  remote  free  stream  velocity  - 
relative 

Air  density  -  lb/cu.  ft. 

Air  viscosity  -  lb/ft.  sec. 

Maximum  value  of  X-inches 
Maximum  value  of  Y  -inches 

Initial  Y-value  for  first  water  droplet -inches 

Vertical  increment  between  water  droplets  at 
initial  X-value  -  inches 

Initial  X  -value  for  water  droplets  -  inches 
Centerbody  highlight  X-value-  inches 
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^  Ysl 

9 

Ysh 

VXWIN 

VYWIN 

DIRECT 

RANDOM 

XHL 

YHL 

ALAMDA 

RES3X 

fv 

Langmr 

GRAVTY 


/I  XI 
AX2 
XREF1 

XREF2 


r 

Minimum  model  Y-value-inches 

Maximum  model  Y-value  for  droplet  calculations- 
inches 

Initial  droplet  velocity  in  the  X-direction- ft/sec. 

Initial  droplet  velocity  in  the  Y-direction-ft/sec. 

Direction  of  droplet  travel  key  -  negative  is  for 
vertical  travel  and  positive  indicates  horizontal 
travel. 

Random  input  key  -  negative  is  for  normal  input 
(increment  Y  drop)  and  positive  indicates  random 
input. 

Cowl  highlight  X -value -inches 

Cowl  highlight  Y-value-inches 

Plow  field  type  key  -  zero  indicates  plane  flow 
and  one  indicates  axially  symmetric  flow. 

Solid  particle -surface  resilience  factor-ratio  of 
the  particle  normal  velocity  after  impingement  to 
the  normal  velocity  before  impingement  (absolute 
values) 

Particle  density  -  Ib/cu.  ft. 

Water  droplet  distribution  summation  key  -mean 
droplet  size  of  Langmuir  D  distribution  if  desired, 
otherwise  input  zero. 

Gravity  and  buoyancy  effects  key  -  zero  if  effects 
neglected  and  one  to  include  effects  in  the 
calculations . 

Coarse  increment  of  X  for  calculations -inches 

Fine  increment  of  X  for  calculations-inches 

X-value  at  change  from  coarse  to  fine  X-increment- 
inches 

X-value  at  change  from  fine  to  coarse  X-increment- 
inches 
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YTRAP1 


Lower  limit  of  particle  trap  -  inches 


T 


YTRAP2 

XTRAP1 

XTRAP2 

TABLE 


Upper  limit  of  particle  trap  -  inches. 


Left  limit  of  particle  trap  -  inches. 
Right  limit  of  particle  trap  -  inches. 


Tables  1,  2  and  9  generation  key  -  calculate 
if  1.0  and  omit  if  zero. 
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Output  of  the  water  droplet  trajectory  computer  program  consists  of 
tabular  print  out  and  a  plotting  magnetic  tape.  The  tabular  print  out 
includes  a  listing  of  the  tables  ana  input  data,  droplet  starting 
conditions,  a  step-by-step  trajectory  record  and  impingement  data. 

If  a  droplet  size  distribution  summation  is  requested,  the  results 
are  printed  at  the  end  of  each  droplet  size  calculation  and  at  the 
completion  of  the  total  distribution.  The  tabular  print  out  includes 
the  following  items. 


XNKL 


ML 


DTW 


CK 

REN 

VEX 


X-value  at  completion  of  calculation  across  an 
interval.  This  value  is  the  XNN  (starting) 
value  of  the  next  interval  -  inches. 

Y-value  at  completion  of  calculation  across 
an  interval.  This  value  is  the  YM  (starting) 
value  of  the  next  interval  -  inches. 

Time  interval  required  for  droplet  to  cross 
interval  from  XNN  to  XNN1.  Units  of  feet  at 
remote  free  stream  velocity. 

Drag  coefficient  parameter,  Cp  Rre^/24 

Relative  Reynolds  Number  of  water  droplet  to 
air,  Rrei. 

Normalized  average  horizontal  air  velocity 
between  XM  and  XML. 


VEXWN1 


Normalized  horizontal  water  droplet  velocity 
at  XNN1. 


VDYWN1 


Normalized  vertical  water  droplet  velocity  at 
XNKL. 


VDY 

XKT 

DIRECT 

WATRAT 

AYS 

XP 


Normalized  average  vertical  air  velocity  between 
XNN  and  XNNl. 

Iteration  of  mean  velocity  counter  used  to 
control  non-converging  solutions. 

Direction  of  droplet  travel  key. 

Rate  of  water  impingement  between  last  two  drops, 
Ib/min-sq.  in. 

Average  surface  distance  from  highlight  for  impinge¬ 
ment  of  last  two  droplets  -  inches. 

X-value  at  droplet  impingement  -  inches. 
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Y- value  at  droplet  impingement  -  inches 


YP 

WATER  Water  contained  in  area  between  starting  points 

of  last  two  droplets  lb/sec.-  in. of  depth. 
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7.0  RESULTS  AND  DISCUSSION 


Comparison  of  the  computer  program  results  with  theoretical  and  NACA 
analytical  and  test  data  for  program  verification  is  shown  on  Figures  1 
through  8.  Figures  1,  2  and  3  verify  the  accuracy  of  the  potential  flow 
computer  program  for  both  plane  and  axially- symmetric  flow  fields. 

Stream  lines  calculated  by  the  updated  potential  flow  computer  program 
around  a  cylinder  at  right  angles  to  a  uniform  stream  (plane  flow)  are 
shown  on  Figure  1.  Theoretical  stream  line  points  are  shown  for  com¬ 
parison,  This  Figure  indicates  complete  agreement  between  theory  and 
the  program  results  in  the  area  that  would  affect  water  droplet  impinge¬ 
ment  with  a  slight  variation  further  from  the  cylinder. 

Figure  2  presents  the  potential  flow  stream  lines  for  a  sphere  in  a 
uniform  stream.  These  data  for  the  axially-symmetric  flow  field  compare 
favorably  with  the  theoretical  stream  line  points  and  are  adequate  for 
accurate  water  droplet  trajectory  calculations. 

The  surfa'  e  velocity  calculated  by  the  potential  flow  program  is  compared 
to  test  results  (reference  22),  on  Figure  3  for  the  front  section  of  an 
NACA  65i-212,  72  inch  chord  airfoil.  These  data  show  good  agreement  in 
the  area  of  water  droplet  impingement.  Disagreement  on  the  aft  section 
is  due  to  intentionally  abbreviated  program  input  and  could  be  Improved 
if  required. 

Figures  4  through  8  verify  the  accuracy  of  the  water  droplet  trajectory 
program  for  both  plane  and  axially-symmetric  flow  fields. 

The  bounding  (tangent)  water  droplet  trajectories  calculated  by  the 
contracted  computer  program  and  from  Figure  2  of  NACA  report  No.  1215 
(reference  2),  are  superimposed  on  the  cylinder-in-a-uniform-stream 
potential  flow  field  on  Figure  4.  These  data  do  not  agree  but  the 
differences  can  probably  be  attributed  to  calculation  differences.  The 
NACA  data  do  not  consider  the  effects  of  gravity  on  the  droplets  and  the 
differential  equations  of  motion  were  linearized  at  distances  greater  than 
five  radii  ahead  of  the  cylinder.  The  contracted  digital  computer  solu¬ 
tion  considers  gravitational  effects  and  is  computed  from  20  radii  ahead 
of  the  object.  The  difference  in  the  bounding  trajectory  ordinate 
values  at  the  origin  results  in  the  computer  program  calculating  ten 
percent  less  impinged  water  than  the  NACA  data  indicate. 

The  water  impingement  efficiency  comparison  on  the  NACA  65q-212  airfoil 
is  shown  on  Figure  5*  These  data  indicate  very  good  agreement  in  curve 
shape  and  near  agreement  in  magnitude  between  the  NACA  TN  3839  (refer¬ 
ence  17),  test  values  and  the  computer  program  results.  A  slight  differ¬ 
ence  in  droplet  size  distribution  from  the  assumed  Langmuir  "D"  distribu¬ 
tion  could  account  for  the  difference  in  the  curve  values. 
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The  computed  water  impingement  results  for  the  l8  inch  diameter  sphere 
are  shown  with  NACA  test  and  calculated  data  on  Figure  6.  The  computer 
program  results  for  the  14.7  micron  volume  mean  droplet  diameter  drops 
compare  with  the  NACA  test  results  as  well  as  the  NACA  calculated  data 
does,  but  they  are  on  opposite  sides.  NACA  TN  4092  (reference  12) 
noted  that  the  accuracy  of  determining  the  test  droplet  size  was  t  6$ 
and  that  this  variation  would  account  for  the  apparent  test  and  calcula¬ 
tion  results  discrepancy.  Refining  the  calculation  increment  (shorter 
time  interval)  to  provide  the  most  accurate  finite  difference  solution 
of  the  differential  equations  of  notion  resulted  in  a  lower  water 
impingement  efficiency.  It  is  believed  that  the  earlier  NACA  calcula¬ 
tions  were  not  refined  as  well  as  the  contracted  computer  program  and 
therefore  were  not  as  accurate. 

The  computed  water  impingement  intensity  on  a  supersonic  engine  inlet 
is  compared  to  the  NACA  test  results  (reference  9)  on  Figures  7  and  8. 

The  centerbody  impingement  intensity  shown  on  Figure  7  indicates 
excellent  agreement  between  the  test  and  calculated  data.  These  results 
show  that  the  computer  program  provides  an  accurate  determination  of 
the  water  impingement  intensity  on  an  axially-symmetric  model  when  the 
proper  inputs  are  supplied. 

Figure  8  illustrates  the  more  difficult  conqiarison  of  test  and  calculated 
impingement  on  the  supersonic  inlet  cowl.  These  data  show  very  good 
agreement  in  both  magnitude  and  location.  The  computed  results  would 
definitely  provide  an  adequate  design  solution. 

Machine  plotting  results  are  illustrated  by  Figure  9  for  the  potential 
flow  program  and  on  Figures  10  and  11  for  the  water  droplet  trajectory 
program. 
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SYMBOLS  AND  REFERENCE 


Symbols 

Symbol 

a 

c 

CD 

D 

F 

g 

k 

K 

r 

» 

R 

t 

U 

U'co 

/V" 

Ar 

X 

y 

/< 

P 

r 


Description 
water  droplet  radius 
constant  of  integration 
Drag  Coefficient 
Drag-  lb. 

Force,  lb. 

Gravitational  acceleration,  32.1?4  ft/see.-sec. 

Proportionality  constant  between  air  velocity 
and  stream  function  that  is  eliminated  by 
arbitrary  stream  function  values e 

inertia  parameter,  2  fa  ^  • 

radius,  in. 

Reynolds  Number 

time,  sec. 

free  stream  velocity,  ft/ sec. 
potential  flow  free  stream  velocity 
local  velocity,  ft/sec. 
normalized  local  velocity,  ZlT/U 
horizontal  distance,  inch 
vertical  distance,  inch 
viscosity  Ib/ft.-sec. 
density  -  lb/cu.  ft. 
stream  function,  relative 
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Subscripts 


a 

i 

ic 

1 

N 

H  +  i 
o 

oc 

ocb 

r 

rel 

w 

x 

y 

©o 


air 

inside 
inside  cowl 
left 

value  at  current  coordinate  location 

value  at  next  coordinate  location 

outside 

outside  cowl 

outside  eenterbody 

right 

relative 

water  or  water  droplet 
in  horizonatl  direction 
in  vertical  (or  radial)  direction 
remote 
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PLATES 


Figure  1 
figure  2 
figure  3 
figure  4 
figure  5 
figure  6 
figure  7 
Figure  8 
Figure  9 
Figure  10 
Figure  11 


Potential  Flow  Field  for  Cylinder 
Potential  Flow  Field  for  Sphere 

NACA  651-212  Airfoil (Surface  Velocity,  4°  Angle  of  Attack 
Water  Droplet  Trajectory  Comparison  on  Cylinder 
NACA  65i~212  Airfoil  Water  Impingement,  4°  Angle  of  Attack 
NACA  TN  4092  Sphere  Water  Impingement 

Supersonic  Nose  Inlet  Centerbody  Water  Impingement  Comparison 
Supersonic  Nose  Inlet  Cowl  Water  Impingement  Comparison 
Potential  Flow  field  for  18  inch  Diameter  Sphere 
Water  Droplet  Trajectories  -  5.77  Micron  Drop 
Water  Droplet  Trajectories  -  18.6  Micron  Drop 
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10.0  APPENDICES 

10 . 1  Appendix  1  -  Potential  Flow  Computer  Program  Flow  Chart 
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3 
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THC  JkZcJ  K/GJT  COM  MW 

AIRf«ANC  Oi  VISION  -  WTCHTTA  tJUXX 


SUBROUTINE  CNSTNT 


srAR.T 


BA-  V(HTrn)+XA 

HME(HON£H)TW) 

BB=  EXrwoiOtic. 

ON£K(£>NE(<+TVVOK) 

pa  c  - 

HTW0(H0H£-t-Hrw0)  ' 
nry.  -F(0NLK.)-ir* 

TY/OKCONEKtT WOK) 

BF  s  — _ t _  •  *' 

D  BAfS&tBCtBDrF 

Yar+4)=-G(BE) 


Y(lT)=Bb*&£ 

1^5 


BA=8A*Y«Nf3)#.|?£ 
Y«T)  =0 

Ycrrt4;=BA+Y(m4) 


J!K, 


rmr+/)»Bs*8£: 


88=88i-YWWt3)¥0£ 
YUTtl) *<?..-  /  -v. 
Y«t-m)  --RB+mr-H)) 


UK 


Yar+2)=BCfBE 


.  \  \  r 


.  (A 


BC=BC*Y«N+0*S£ 

YUT+D--0 

YQT-f 4)=BC  +YCIT-f-4) 


E  1666  RS 


HCVUtO 


P 
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*  vC“/»  *  4‘**»  '  •  •  ’ 


► 


I 

t 

I 

I 

{ 


CNSTNT  X 


™  £?g?j?s/y&  < 


AmPLANC  OfVlftlCN  -  WlCWTTA  »*AJ* CH 


*  *v- 
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-MN 


SrO/?£CJH)»Y(WLOC)] 


ks 


£$P  4* = (RCLAXCNLOCJj  NLCC } 
~Y(NLOC))*0Ni££A 


'write,  mjc 
values  stkt7n$ 
wth  ycncfct) 

ON  Taps  t  AND 


vmzwjji: jj 


write¥swee?  num$£.r  my  residue. 

SEQUENCE.  (RSPl 

COLU'/.N  {HFU)j  ROW  (NFUN)}  WITH 
CONVERGENCE  HISTORY  (STGZCCnW.J.N 


0/  ‘  a  -*•  , 


inX^sv;^ 


NJW^:  Hi 

l^alHi! 

ZNSWPS 


NSWE.P1  zNSv/UP 
N*Z  ' 


JNl*JNltl 


H 


V  • 

V*  '  i  - 


S\V££P,  4 


*■■-** •***’*-*'  =t  t-  ff*.y  ^ 


OfltlASM 


THE  COM  PANT 

AI»PLANC  01  VISION  *  WICHITA  BKANCH 


r  i 


•=) 


i  * 
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3 


NBE7AD=NBETAC-1 
NNYR=NCFCT  • 

II=K  J 

’ 

V 

4. _ 

V 

t _ 

NNYR1=NNYR+NYR-1 

NNYRc=NNYR+t 

NNYR3=NNYR1+1 

;;  . :r 

/^BETA^\  >0  XNNYR-NNYR+I  [  J 

WRITE  "PROGRAM  LEFT 
SUBROUTINE  SWEEP  AFTER 
(IN)  SWEEPS  WITH  <RSD*M 
RESIDUE** _ 


I=NGRDXC] 


) 


3 


frci  >=o 

i=x+i_ 


firm  ovc 


<t> \  .. 


/ 


(RETURN^ 


ftHt  I.Z 

e  1606  ns 


SviESP  6 

saBsgaaBBsssEBagEGsaasBaaBaBggaga^^^^ 


» 
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FUNCTION  RELAX 


^  •  • 


2 


fr»  *v?  •JfffA '  .  .-•jn’W'N 


t  1663  Ri 


nzaesagasasgc  jaagsHfflBsaBSBBSs^sEEaBBSBB&^SBKSu  ^gasssso. 


SS355SK®J 


i 


t 
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THE  COMMN Y 

AMnjtNt  DIVISION  •  WICHITA  IMNCM 


E  1 686  Hi 


n, 


« j&&jsf£mep  COMWW 

AI«n_ANC  CSV* •  ION  •  W1CHCTA  SlUNCH 


D3'  6961 
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NfUNCr*NCFCf 

NmAPx^mRoxc 
NMAP'fzUGRPre 
N6ZAD'S  =  MQR£M: 

u — 7ZEZT - 

KP*l 

/wRlTE.  *X*  \ 

J  (YCJ-)?"anO  j 

l  COLUMN  J 
\  ffgAPMGS  / 

W57>/ 

KPX*/ 


/Wir£  ya), 
Y(nfunct)  Y 
WG£40tyy(Nj*Al 
^Y(N6RAt7V)  * 


'wRjrc  ycd 

2.£RO,  iwo/ 

ZCRO,  2£R0 


HFUNCT-NFUNCr+i 

NSRApy-N^RADYtl 

NGRADX-NGRADX+t 

NGRADV-NGRAPVfl- 

m*KPx+i  . 

«5>^ 

pNYR 

JWNT(Y(J+I))  +Z 
KP-'KPH  _ 

<i>^ 

|>NXC 

^  RETURN^ 


V^ZQC  X 


■Witenar 
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COM  WW 


orvi»iOH  -  wichtta  *junch 


SUBROUTINE  ICMP6S 

«  *  *  . 

START  *• 

_ : _ * _ : _ 

W6RADV-NSKDVC-/ 

NCf  --  .. 

NXYAA-I  c 
N6RADX=N6RMCi 
PSPADY.-HGIWC-I 
UKXVtNKXVt-l'-  • 
NFUNCT-NCrcT-l  . 
NDX*NPXC-J. 

S-KJ' 

®lii£rj  J  . 

T:  — | 

ND^NDYC-/  .  .• 

N  K.YV=NkYUCd  ‘ 
i-  NOX-NO X+|j-.  . 

IH5T+I.'  •• 

.  -  MKXy  =  NKXV+f 

r.j_L, 

„  .  |nY»(  I  , 


NGWX-NSRfiPX-N; 

NtfRADY  =NGRAPY+/ 
N6KA0V-’NGRA0\'+I 
N  FVNCT ~WFUNCf+l 
NKYY-NKYVvf, 
Ndv-NDY+1 

p;xz, 

tTN-J+a. 

IN  *1+2' 
KYVAsYCNKVV>4 
KYVG--  KYVA*  ‘ 


RCVXILO 


T)4C  £3i COMPAW 

AI*»VA>«C  Oft.'»W«»  «  V/lCIUTA  VftkMOe 
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SUBROUTINE  GRDSTI 

:  M*£L 
, — ~I_ — , 
sMXM- 

Mi  A*! 

IN  f  -NFUNCP f  HYR 

|N33  NFUNCTft  nyk  .  .  '  a  ' 

Hf*Y(NPXy!  ‘  ' 

H2*Y(MPX-l)  '* 

77Ci=Y(W0Y)  .  . -  • 

TKX*Y(WJiH)’ '  .  *  '•  * 

UO=Y(NFUNCT) 

ut^i(m)  *  ■  . 

UZ=Y(NrVNCT+l) 

03-YCN3) 

U4  =Y(N  FUNCTr# 

I<T/<=0  •  J'  ..  .  ♦"  • 

ikj*0 

KIJc0  :  \  •.  ;  • '  . 

NIK»0  1 


'Y(IN)'Y(fr\  >0  v  tjJWO  =Y(ThiaN)  - 

sYcndx-;)/  7 U3*Y (IM+O . .  • 

\  x  IJK-I 


YCIN-+2.} 

~YCJ)*YM 


HQNE:*Y(WtZ)-Y<t 
VI  *YON+$ 


YM-YU> 

+YCN0Y4 


meKzya)~Y(jN) 

u$~-rmtO  •  I 

KIJ"  *  I  ! 


RETURN 


-YdhYfml 


oMEK'YUNtzynv 

VZ.-Y(0Nt3). 

JIM _ _ 


TWt  GMJSTMJ  \SG*r  COMPACT 

AINP1UNC  DIVISION  •  WtCKCTTA  BKANCX 
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SUBROUTINE  SUBB1 


r  '>  nt,  i  «*  .**•»  :* v r»  >>r^\  *  /♦.rs*  t  f  >v  ?tV~»KiT4  *  I yy»»*>  ***  v  .*•  y* 


ncvuxo 


E  1686  AS 


COM  MW 

AinruANt  Division  -  wtcxtta  maxh 


THC 


COMMW 


airpunk  division  -  wicmita  branch 
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© 


NPWKAD'f-NW 

NXA*J 

NJ  -NFUkJCT-NYK 
W3^MI-N  yR' 
U0--YC/Nf3) 

Ul  *  Y(N/; 


\m=YCNPx-i) 


t;3*YaNri)  ~ 


U3-YCH3) 

H2’Y(NDXrX) 


€. 


\h 

\ 


— V*  ‘  • '  -  H. 


V  '  / 


rr-Y(MDr-i) 

7%  -  YCNOYi*/; 

uz- Y(nfunct-h) 
Ub*Y(NFUNCTt3) 
U4-Y(NFUNCT~1) 
us'-Y(NFwcr-i) 
T*X=Y(NDY) 
tkx-Ywy-i) 


t  1685  R» 


I 


CO 


O  s 


N3 = NFIWnI  CT  -M  VR 

nya=3 

V0-Y(iTN+3) 

Vl'Y(NFUNcT-l) 

TKJ-Y(NDy-l) 


rf6myW^M75jU# 


U4=Y(JNH)  ■  - 
TKl  =  |Y«)-Y(JN)j-TXl 


ncvtSKO 


0 


com  mw 

AIAPt_4NE  OlVietOH  .  WICHITA  *HMCM 


1  >  J 
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SUBROUTINE  DERJV 

STUART  . 

□r  ' 

(fcT=rj  ■ 

<^yip^zL ••- 

?/  Fct=Y03 

^ - - — . i 

El  ED 

- 2 

paTA»)l  • 


06LTA=| 


HZ=tfH-H2. 


*  / 

Vk  / 


,  DELTA (H2.2- Cui  ~U0)-h(HI2’{uQ-u3) ) 


A  (HI)  (HZ)  (H 2-KFS)  (HD) 


=  S 


RETURN 


-GX 

VaV 

*1 

f«4 

e  t{««  r* 


COMPiW 


AIMPUNC  OIV19ION  -  WICWTTA  tKAMCM 


SUBROUTINE  LINES 


D  3 

Page 


JLN'ILN+I 


UN  =0 
tiNL  NS*I 


RETURN 


READ  NSlj  y 

cscrij  With, 

V 0  i 


c$  CH 


N6RPC-NSRDXC 


NNLfc 


lwNLNS~l 


(REAP  ’  \ 

sunescjnk;,] 
JN^I,N5l-  J 

h^Z'ZCIZI, 

NSLNC=A/GRpC 

zzxi — 

NS*  I 

r—E^r 

LN-MSLtfC 

nizL—J 

(  MARC  \ 


N6RA0-NGPJX. 


IflfeHgfegi 


3  ') 


3) 


mxvicko 


acviftfLO 


ccMiwr 

AIUPUM  DIVISION  -  W1CKJTA  CftAXCH 
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WCV  l««LCl 


COMPANY 

AIKPUNt  OtMtlON  •  WICHITA  IKIXH 
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A  . 


NSTOK-NSTOR 

NST0R-LN‘+6 


UN£*0 

XUP-XLl 

yup^ix 

YUP=YU 

YL2p*YL2. 


1 

1 

V 

VP-  l'CIt3l 

'  '  VO 

write  Y(i-H).rat2),  ^ 
Y(it3),  vA,rc  M?,Y(w>y 


THE 


J&4P&SA Y& 


COMMNT 


AIRPLANE  DIVISION  •  WlCHfTA  ■  RANCH 
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m  cvtuo 


.)■ 


ncvtsKo 


THE  \S iU-  COMPANY 

AIRRuANC  OIVIS-ON  -  WICHITA  •  RANCH 
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e  166tr  Rg 


lAdBU  ,  u  aAv 


*CVl*EO 


T  1  -j^ry 


2S5SSE5SZSS 


THE 


ME? 


COMMNY 


*mrj>Nt  Division  •  wichita  iiunch 


A 
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SUBROUTINE  MARC 


hl-0 

smK-i 

FLA6=0 

2=0 

LU'NSRDC 

i\'6.RADXr  ti$R  DXC'  I 

NGRADY-N6RDYC-J 

A/V-NSRDVC-I 

NKXV^MfcXVC-/ 

MFUNCT^NCFCr-l 

npx^hdxc-i 

J'KJ 
FST-0 


MX*  I 


MA/K= SPANK 
l!Dy-Noyc~J 
HKYV^^KyVC-l 
UDX-UpXfJ  . 

i-isr+i 

NW-Nk/jV-hl 


ILN-'lN 

N6RA0X=M(?/^P/+/ 

NG^^Y-MG.^PY-f) 

N'RUNCT^MFMCT-H 

NV-’NV-H 

NKYV-'NKY  U+/ 

UPY'NPY-tl 

JN-J+X 

kywyowvw 

KYV3=  KYM 


E  1686  at 


3Z55E 


QiflASU 


«cvisro 
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HCvuro 


Hcvuro 


™  Z7^7Z50Mfy  COM?W 

OtVISlON  -  WICHITA  5^A>VSW  ' 
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|  kpnt*ln+(,  -ik 

j  YC<Otir)  “Y (KGNT-s) 

i  K>rx+i _ 


W~ 

T>1XK 

Yttox)-za 


^  IKK 


2.5  1 


^  X 

R4.VO - 


5£/\NK-2.,0 

RANK»SMNfc 

fST=  /.(? 


K.ONT-  LMf  6-1 K 
YCKG,vd-Y(K0/lT-y; 
l£zlK±l _ 

Y(ui-w)-i 

raA/-n>3 

Y  (IN- 20)- 1 
X(LM-lf)*4 
Y(laj-io)--l 
YCLH-f)-2> 
YUK)  -  i 
LN-LA|4!T _ 


KCNT^LNi'b-lK 

YiKoNT‘)^YOCJjHT-5) 
IK -I Kt/ _ 


YCLN-i5T)»Z 

YON  ~20)  O 

YCuNS) 

Y(jJi'\o)-x 

yQjJ'*r);J> 

Y(<-H)s1- 

iM-LbifS' 


If 


ii  thc  IZ?£PIjFj7f^f/[iT , 


:—r.  ^ 

«  I  ■  ■■■  ■  «  ».»  -  *  *«■ 


COMftVN/ 

AIUPUANC  OIVtSIOM  -  WICHITA  SRANCM 
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Appendix  2  -  Potential  glow  Computer  Program  Listing 


NO.  D3-6961 

REVLTR: 

SECT 

PAGE  -112 

K  *30  S  8  R» 


O  */* 


I8FTC  W AMAIN 

MAIN  PROGRAM 

DIMENSION  SiU.YC  18000) , SL INESI 30) 

COMMON  6,NXC*NYR*NKXVC,NKYVC»NDXC»NDYC  »N'Cf  CT  *N8ETACjKJ»8I TS* 1ST, 

I  HONE  ,  HTWO ,  ONER ,  TWOK ,  I »  !  N,  J ,  JN  *  I K  J ,  K I J ,  J I K  *  1  JSC ,  i  T  ,  NDX  ,  NDY  *  NFUNCT  ,  Y,  0* 

2SL I NES, HI »H2, TK1 » TK2  »N8 1 »NB2  »N83  *N84»NC l »NC2 , NC3  »NC4  »Ni »N2*U0,Ul »  0( 

3U3»U2,U4,  1PRIN»NRS0UE»NSL, NSLNC,NGRDC ,NGRDXC  »NGROYC, AXYM,NGRAO*  (H 

4NGRACX  , NGRAQY, LN,CSCH,NG»NDER I V, SSRCH, NSWPS, OMEGA, CVGS, CONST,  0( 

5STREAM,KNTR, A, C, 0, E , F, G» HK i tHK12 , H12, HK2 , 3 Y,BX , NKXV, NKYV, N3,  0( 

6NLRCAL  »Q,RfS»T»V,W»XC,YC»IH,NTAPE  OC 


COMMON  /VEL/  NGROVC 

COMMON  /C0M2/  XSCLtYSCL, PUNCH,  I TL6 , Y0? YI , XL l , XL2 , YU * YL2» XSCIA, 
1YSCLA»XL1A,YL1A,XL2A»YL2A 
DIMENSION  I TLE ( 24 ) 

REWINO  23 
l  DO  10  1*1,18000 
10  Y( I ) *0# 

S  =  1 » 


PUNT  *  0,  0( 

HEDGE*. 75  0( 

NSWEEP*20  0( 

FL I M= 10.  0( 

CALL  INPUT 

902  CALL  SEARCH! Y0,Y1,X,KI)  0( 

C  THIS  ROUTINE  ALSO  FINDS  ERRORS  IN  THE  INPUT  DATA , GENERATES  THE  0( 

C  COEFFICIENTS  OF  THE  DIFFERENCE  ANALOGUE,  AND  ENFOFCES  DERIVATIVES.  0( 

IF(KNTR)40,907,40 

907  CALL  SWEEP(PUNT,NSWEEP, HEOGE , FL l M )  Q( 


V*I. 

CALL  SEARCH ( YO, Y l » X »K i ) 
CALL  V.ELOC 
V  -- 1 . 0 

CALL  SEARCH  (Y0,Y1,X,K1) 
CALL  LINES 
40  STOP 
END 
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itBFTC  WALNKi 

SUBROUTINE  INPUT 

D I  PENS  LON  Y  (  18000} «  SLINESC  30)  ,X(6) 

COMMON  S,NXC,.NYR!.NKXVC,NKYVC,NDXC,NOYC,NCFCT,NBETAC»KJ,BITS,IST, 
1H0NE  ,HTWO»ONEK»  TWOK  ,  I  *1 N ,  J ,  JN ,  IK  J ,K  I J,  J I K *1  JK ,.1  T  , NOX , NOY , NFUNCT,  Y* 
2SL1NES,H1,H2,-TKI,.TKZ,NBI,NB2,NB3,NB4,NC1,NC2,NC3,NC4,NI,N2,U0,UI, 
3U3,U2»U4,IPRINj.NRSDUE»NSL.NSLNC»NGRDC,  NGROXC  »NGROYC , AXYM*  NGRAD, 
4NGRA0X»NGRA0Y,LN,.CSCH»NG»NDERI  V,  SSRCH,NSWPS  » OMEGA »CVGS» CONST* 
5STREAM,  KNT  R,A,C,0»E,F»GtHK>i,HKi2».H12,HK2»BY  »BX»NK.XV»NKYV»N3  » 
6NLRCAL»Q»R»S,T».V*W»XC»YC»IH»NTAPE 
COMMON  /C0M2/  XSCL,YSCL, PUNCH,  I  TIE  » Y0» Yl *  XL 1 *XL2 »YLi »YL2 » XSCLA, 
1YSCLA,XLIA,YL1A,XL2A,YL2A 
COMMON  /V£L/  NGRCVC 
DIMENSION  ITLE124) 

501  FORMAT  {6E12.6t 

502  FORMAT  (6112) 

503  FORMAT  (12A61 

600  FORMAT  (7H  AXYM  =F2 . 0,.!  IX,  7HNDER IV*  1 2‘,  10X,5HXSCL*F8.3,4X,5HYSCL*. 

1F8.3/  7H  NSWPS=I6,  7X,.7HY0  =F8.3,  4X,5HXL1  *F8 .3 ,4X,5HYL1  * 

2F8 • 3/  7H  CVGS  =F9.7,  4X,7HY1  *F8.3,  4X,5HXL2  *F8*3,4X,5HYL2  * 

3F8. 3/  7H  PUNCH  =  F2.0,30X,6HXSCLA-F7.3»4X,.6HYSCLA*F7.3/ 

4  39X,5HXL1A=F8.3^4X,5HYL1A=F8.?/  39X5HXL2A=F8.3 ,4X» 

55HYL2A=F8.3//) 

610  FORMAT!  1HU12A6/1X,  12A6///1 
XSCLA*0. 

2  A*=l. 

V  —0* 

CSCH=0. 

N2*0 

Tfc0, 

C*l. 

CONSTfc l# 

CM6GA- 1 . 37 
STREAM  *  o„ 

RUN  -  0. 

DUMI*0« 

READ  (5,503)  ITL6 

WRITE  (6,610)  ITL 6 

READ  (5,501)  AXYM, CVGS,ANSWPS, PUNCH 

NSWPS* ANSWPS 

READ  (5,501)  AND6RV , YO*Y i 
NDER l V=ANDERV 

READ  (5,501)  XSCL,YSCL,XLl,XL2,YLi,YL2 
READ  (5,501)  XSCLA  ►YSCL  A  ,XL1A,XL2A,YLIA,YL2A 
IF  (AXYM, EQ. 1.1  STREAM^- U 
•IF  (  AXYM, EQ .  I •  K  C0N$T*0. 

900  CALL  CHECK 

IF  (KNTR.NE.O)  STOP 

WRITE  (6,600).  AXYM,  NDER I V ,  XSCL,YSCL,NSWPS,.Y0,XLI,YL1,  CVGS ,  Yl  ».XL2  , 
1YL2, PUNCH, XSCLA, YSCLA,  XL1A, YL 1 A,  XL2A, YL2A 

42  RETURN 


Of 

0( 

0( 

01 

0(, 


0(1 


0( 

0( 

0( 

o(: 

I 

0( 

0(, 

0( 


! 


0( 
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D3-69& 
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o  o  o 


SIBFTC  WACHCK  \  . 

SUBROUTINE  CHECK 

C  THIS  ROUTINE  READS  COORDINATE  INPUT,  CHECKS  FOR  ERRORS, AND  PRINTS 
C  COORDINATE  INPUT. 

DIMENSION  B( 1),Y< 18000) »SLINES(30) 


COMMON 
1H0N5 i H 
2SLINES 
3U3 » U2» 
4NGRADX 
5STREAM 
6P,Q,R, 


>SfREAM,KNTR, A,C, D, E , F, G, HK 1, HK 12, HI 2»HK2 , BY,8X,NKXV, NKYV,N3* 
>P,Q,R,S,T,V,W,XC,YC»IH,NTAPE 
COMMON  /VEL/  NGROVC 

FORMAT  (4H  Yl=F12.5/6lH  CHECK  THE  INPUT  DATA  CARD  WHICH  HAS  THIS  Y 
(  AS  ITS  FIRST  TERH/47H  EITHER  THE  CONSTANT  TERM  IS  NOT  DIVISIBLE  B 
(Y  4/55H  EXACTLY  OR  THE  BOUNDARY  COORDINATES  ARE  NOT  MONOTONIC.//) 
FORMAT  { 4H  Y1-F12.5/4H  Y2*F12.5/41H  THE  COORDINATE  VALUES  ARE  NOT 
rMnNnmwf  r .  //  » 


ecu  r  urvnm  i  til  i  i  —  r  l  c «  Jf  *tn  i  nc  ouum/IIim  i  C  tmuuco  Mf\C  itUI 

XMONOTONIC.//) 

501  FORMAT  (6E12.6) 

599  FORMAT  f 1HI 1 

600  FORMAT  (67H  COORDINATES  AND  FUNCTION  VALUES  FOR  BOUNDARY  POINTS  ON 
1  MESH  ROWS  //8X, 1HY,13X,1HN,13X,1HX,12X,6HF(X,Y) ,10X, lHX,llX,6Hr( 
2X , Y )//  ) 

601  FORMAT  I 68HIC00RDINATES  AND  FUNCTION  VALUES  FOR  BOUNDARY  POINTS  ON 
l  MESH  C0LUMNS//8X, IHX, 13X, IHN, 13X, 1HY, 12X,6HF(X,Y) ,LOX, 1HY, UX,6HF 
ZlXjY)//) 

602  FORMAT  ( IX, 6F14. 7,/ { 29X, 4F 14.7) )  ' 


L  riLJIl  VULUnH  /  OAf  IflAf  i  jAf  Jl  nil  f  LJ 

2(X'/Y)//) 

602  FORMAT  ( IX, 6F14. 7,/ { 29X, 4F 14.7) ) 
620  FORMAT  (iX,6F14.7) 

686  FORMAT  ti.v>,4IlO) 


v  VA.  '  W  f  >  •  I  »»  • 

620  FORMAT  (iX,6F14.7) 

686  FORMAT  (l.%»4110) 

LINE=*5 
KNTR  *  0 
I STs600 
I  *  ISTfl 

Y { I  )  IS  THE  LOCATION  OF 
NKYV=l 


4  4 

I  *  ISTfl 

YII)  IS  THE  LOCATION  OF  THE  FIRST  Y  INPUT  VALUE. 

NKYV=l 

NKYVC»NKYV. 

NKYV  AND  NKXV  GIVE  THE  BEGINNING  OF  KYV  AND  KXV 
IN  Y  ARRAY.  THESE  TWO  VECTORS  TELL  HOW  MANY  BOUNDARY 
VALUES  ARE  CONNECTED  WITH  EACH  Y  AND  X  COORDINATE. 
NDY= 150 


NDY= 150 
NKXV=300 
NDX=450 
KNT  3 1 

WRITE  (6,600) 

I£«I*5 

READ  (5,501) IY( IR) , IR=f,IEl 
L  =  YII+1) 

Y (NKYV) «L 
40  XL=L 

IF  (AM0D(XL,4.))240,44,240 
240  WRITE  (6,227)  Y(I) 

KNTR  *  i 
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46  LN  *  L-2 

C  LOOP  5  CHECKS  MONOTONICITY  Of  X  BOUNDARY  COORDINATES 
DO  5  M  *  l,LN,2 
IF  (M.EQ.l)  GO  TO  52 
XM*M+l 

IF  (AM0D(XM,4.)  .NE.O.)  G"  TO  52 
I 

IE* IB+3 

READ  (5,501)  (YUR),  IR*  IB,  IE) 

52  IN= J+M-i 

I T { Y  C tN+4)-Y( IN+2))245,5,5 
245  WRITE  (6,227)  Y C J I 
KNTR  *»  l 
5  CONTINUE 

LQ*J  +  INT(Y(J-*-i)  )♦! 

WRITE  (6,602) ( Y( L4) ,L4*J,LQ) 

L  INE*L INE  +  I 

IF  (LINE. LT. 25)  GO  TO  62 
LINE*0 

WRITE  (6,601) 

62  MIT*J+L*2 
IE*M!T+5 

READ  (5,501)  (Y(IR),IR*MIT,IE) 

IF  (YIMIT  +  I )  .EQ.O. )  GO  TO  285 
KNT =KNT-f  l 

Y(NDX)*Y(MIT)-Y( J) 

I F ( Y ( NOX ) ) 243,8, 8 

C  IF  STATEMENT  CHECKS  MONOTONICITY  OF  X  COORDINATES 

243  WRITE  (6,228)  Y(J),Y(MIT) 

KNTR  *  I 

8  NKXV  *  NKXV  ♦  l 
J=MIT 
L*Y( J+l) 

NOX  =  NOX  *  l 
26  Y (NKXV ) *L 
GO  TO  260 
285  NXC*KNT 

00  290  I Q= 1 , NYR 

I PT=NYR* IQ 

Y(  I PT  )  =Y (  IQ-H49) 

290  CONTINUE 

DO  291  IQ*l,NXC 
I PT=2»NYR+ IQ 
Y (  I  PT ) =Y ( IQ  +  299 ) 

IPT  =  2«NYR>NXC4-IQ 
Y (  IPT  )  *Y ( I Q+449) 

291  CONTINUE 
1$T=2»(NXC*NYR) 

NOY=NYRU 

NDYC=NDY 

C  NDYC  GIVES  LOCATION  OF  FIRST  Y  MESH  VALUE  IN  Y  ARRAY. 
NKXV*NDY*NYR 


I 
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44  LN  *  L-2 

C  LOOP  16  CHECKS  MONOTONICITY  OF  Y  BOUNDARY  COORDINATES 
DO  16  M=l,LN»2 
IF  (M.EQ.l)  GO  TO  51 
XM=M+l 

IF  (AM0D(XM,4.)  .NE.O.)  GO  TO  51 
IB*I+H+3 
IE*  18+3 

READ  15,501)  (Y( IR ) , IR=I8, IE) 

51  IN  =  i+K-i 

I F { Y (  1N+45-Y (  IN  +  2)  )241, 16,16 

241  WRITE  16,227)  YII) 

KNTR  *  l 

16  CONTINUE 

LQ=I+INT(Y(  1  +  1)  )  +  l 
WRITE  (6»602)(Y(L4)*L45SI»L<3) 

LIN£aLlN£+l 

IF  ( L INE  »LT«  25 )  GO  TO  61 
UNE*0 

WRITE  (6,599) 

WRITE  (6,600) 

61  MIT  *  1+1*2 
IE*MIT+5 

READ  (5,501)  (YIIR),  IR=*MIT,  IE) 

Y(NDY)*Y(MIT)-Y{  I ) 

IF  (Y(MIT+l).EQ.O.)  GO  TO  280 
KNT  =KNT  + 1 

I F ( Y { NDY ) )2 44, 24, 24 

C  I F  *  STATEMENT  CHECKS  M0N0T3NICI TY  OF  Y  COORDINATES 
244  WRITE  (6,228)  Y ( I ) , Y ( Ml T ) 

KNTR  *  l 
24  NKYVaNKYV+l 
I  =M  S  T 
L  =  Y  1 1  *1 ) 

NDY  a  NDY  +  1 
4  Y  l  NKYV  )=*L 
GO  TO  40 
280  NYRaKNT 
KNTal 

J«2vL+I 
K  J*  J 

C  YtKJ)  IS  THE  LOCATION  OF  THE  FIRST  X  INPUT  VALUE. 

L INE*0 
IE=*J+5 

READ  (5,501)  { Y ( IR),IR*J,IE) 

L*Y ( J+l )  0 

Y ( NKXV ) 3L  0 

WRITE  (6,601) 

260  XL*t 

IF  ( AMODfXL ,4. ) ) 242,46, 242 

242  WRITE  (6,227)  YtJ) 

KNTR  *  1  0 
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NKXVC*NKXV  0 

NDX=NKXV*NXC  0 

NDXC=NDX  0 

C  NOXC  GIVES  LOCATION  OF  FIRST  X  MESH  VALUE  IN  Y  ARRAY,  0 

00  295  1*1,4000 
I9*i$r*i 
YIID)  =*Y  1 1*600) 

295  CONTINUE 

K J*KJ-600*I ST 
NCFCT =MI T-60Q+I ST 

C  Y  I NCFCT }  IS  THE  FIRST  FUNCTION  VALUE  LOCATION. 

NBETAC*NCFCT+(NXC*NYR> 

C  Y(NSETAC) IS  THE  FIRST  CONSTANT  VALUE  LOCATION. 

MK*NXC*NYR 
NGR0VC*N8ETAC 
NGROXC*NBETAC*MK 

NGR0YC*NGRDXC+8K  .  . 

NGRDC*NGRDYC+MK 
IF  (LINE. GT. 22)  WRITE  16,599) 

59  RETURN 

END  0 
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$  1 8FTC  WASRCH 

SUBROUTINE  SEARCH! YO»Yl *  X»K1)  0 

THIS  ROUTINE  PLACES  THE  INPUT  BOUNDARY  VALUES  IN  THE  PSI  ARRAY 
WHICH  BEGINS  WITH  Y ( NFUNCT ) «  IT  CHECKS  FOR  INCONSISTENCY  IN  INPUT 
AND  CALLS  REGION  FOR  EACH  MESH  POINT  INSIDE  THE  REGION. 

DIMENSION  BID,  Y{  18000),  SLINESI 30) 

DIMENSION  X(6) 

COMMON  B,NXC,NYR,NKXVC,NKYVC,NDXC,NDYC,NCFCT,NB£TAC,KJ,BITS,IST, 

1H0NE , HTWO, ONEK, TWOK, I*IN,J,JN, IXJ,KIJ,JIK,IJK,IT  »NDXt NDY , NFUNCT » Y, 
2SLlNES,Hl,H2,TKl,TK2,NBl,N82,NB3,NB4,NCl,NC2,NC3,NC4,Nl,N2,U0,Ut, 

3U3»U2 ,U4, IPRIN»NRSDUE,NSL»NSLNC,NGRDC»NGRDXC»NGRDYC, AXYM»NGRAD* 

4NGRADX, NGRADY,LN,CSCH,NG, I XTRAP,SSRCH,NSWPS, OMEGA, CVGS, CONST, 
5STREAM,KNTR, A, C, D, E , F,G, HK1 , HK 12 ,HI2,HK2 , BY, 8X ,NKXV, NKYV*N3 t 
6P,Q,R,S,T,V,W,XC,YC, IH,NTAPE 
COMMON  UBDY,L8DY 

225  FORMAT  (4H  Y1*F12.6/4H  Y2*F12.6/56H  THE  POINT  (Yl,Y2)  HAS  ONE  COQR 
XDINATE  OUTSIOE  THE  REGION//) 

620  FORMAT  !IX,6F14.7) 

621  FORMAT  (3124) 

625  FORMAT  (120H  THIS  PROBLEM  IS  TOO  BIG.  ELIMINATE  AS  MANY  OF  THE  MES 
IH  POINTS  FOR  WHICH  PSI  VALUES  ARE  CALCULATED  AS  THERE  ARE  OF  SUCH- 
2/26H  POINTS  TO  THE  RIGHT  OF  X*F14.7| 

626  FORMAT  (8F12.5) 

630  FORMAT  (5F15.4) 

631  FORMAT  (IHU 
IF  (CONST. EQ.l. )G0  TO  8 
A*l  • 

C*l. 

0*0. 

G*0. 

F*0. 

8  BITS=t-2.**l5 
UBDY=0. 

L8DY-0 
KNTR  =  0 
I T=NBETAC 
NKXV*NKX VC” l 
NFUNCT-NCFCT-1 
NDX*NDXC”  1 
J*KJ 

C  LOOPS  75  AND  76  PERMIT  CHECKING  POINTS  BY  COLUMN 
00  76  NX* i, NXC 
NDY=NDYC-1 
NKYV=NKYVC”l 
NKXV3NKXV+l 

9  NDX=NDX*l 
XC*Y{ J) 

i  *  ism 

DO  755  NY*1,NYR 
N0Y*NDY+1 
NKYV*NKYV+l 
NFUNCT*NFUNCTU 
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YC*Y { I )  0 

JN*J+2  0 

lN*I+2  0 

KYVA=Y (NKYVJ/4.  0 

KYV8  *  KYVA  0 

LOOPS  10  AND  112  CHECK  ALL  POINTS  BETWEEN  BOUNDARY  COORDINATES  0 

TO  DETERMINE  THEIR  LOCATION  IN  THE  FIELD  0 

DO  10  I  Sa  1 »  KYVA  0 

ENTRY  PERMITS  TRANFER  OF  BOUNDARY  MESH  POINTS  TO  FUNCTION  ARRAY  0 

ENTRY=Y{ IN+il  0 

I F  ( Y  ( IN )  —  Y {  J  ) )  1  >  2»  3  0 

Y ( IN)  IS  THE  LOCATION  OF  THE  LEFT  X  BOUNDARY  COORDINATE  0 

1  2NTRY*Y { IN+3 )  0 

t  ^  { Y  {  lN+?)-Y(J)  )U0, 2, 7  0 

7  KXVA=Y(NKXV)/4.  0 

KXV8=KXVA  0 

DO  112  JS*1,KXVA  0 

ENTRY*Y{  JN-H)  0 

IF(Y{JNJ-YU  mi, 2,33  0 

11  ENTRY*Y ( JN+3 I  0 

I F ( Y ( JN+2 )-Y ( I )  )  12*  2 » 13  O' 

SUBROUTINE  REGION  GENERATES  MESH  DISTANCES  THAT  0 

ARE  NOT  FOUND  IN  THE  CHECK  SUBROUTINE.  THESE  ARE  0 

THE  SHORT  DISTANCES  TO  THE  BOUNDARY.  0 

13  IF.  (V)  75,113,75  0* 

113  CALL  REGION  O1 

GO  TO  75  Ol 

12  KXVB=KXV8-1  ’  0 

IF(KXVB)  33,33,112  0 

112  JN= JN+4  O1 

110  KYVB  *  KYV8-1  0 

IF  (KYVB)  3,3,10  0 

10  I N* I N+4  0 

Y(JN)  IS  THE  LOCATION  OF  THE  LEFT  Y  BOUNDARY  COORDINATE  0 

IF  A  POINT  IS  IN  THE  FIELD  IN  THE  X-DIRECTION  0 

AND  OUT  IN  THE  Y-DI RECT ION,  THE  ROUTINE  GOES  0 

TO  LOCATION  750  INDICATING  INCORRECT  INPUT  DATA.  0 

3  IF  (Y(JN)-Y(I)) 101,750,333  0 

101  IF  (Y(JN+2)-Y(I ) )333,750,750  0 

750  WRITE  (6,225)  Y ( I ) ,Y{ J) 

KNTR  *  1 
GO  TO  75 

C  MINUS  ALL  BITS  ARE  STORED  IN  COEFFICIENT  ARRAY  FOR 
C  POINTS  WHICH  ARE  NOT  IN  THE  REGION  DEFINED. 

33  IF  (Y{ IN)-Y(J) ) 171,750,333 
171  IF  ( V C IN+2)-Y(J) ) 333  ,750,750 
333  IF  (V)  400,633,632 
633  Y { I T ) =B  I T$ 

I T=IT+l 
GO  TO  75 

400  IF  (Y(  JNI.GT.YU  ) )  Y  ( NFUNC  T )  *Y(  JN+1 ) 

IF  ( Y { JN+2 ) • LT. Y ( I } )  Y ( NFUNC T ) *Y { JN+3 ) 
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GO  TO  75 

632  Y (NFUNCT)*BIT$ 

GO  TO  75 

2  IF  (V)  75,22,75 
22  ASSIGN  2222  TO  Ll 

IF  ( IXTRAP)222,232,222 
232  CONTINUE 
2202  Y I NFUNCT ) *ENTRY 

GO  TO  Lit  12222*75) 

2222  YdTHBITS 
IT  *ITH 
GO  TO  75 

C  CHECK  HERE  TO  FIND  WHERE  DERIVATIVES  ARE  TO  BE  ENFORCED 
222  IF  (Y{ I1-Y0J200,  15»200 
15  IF  (YtJJ-Yt INJJ25, 232,25 
25  IF  { Y  C J)-Y{ IN+2) )123,232,123 

200  IF  (Y(I)-Y1)232, 321,232 
321  KN*1 

DO  201  K3 1 » K1 

IF  (Y(J)-X(KN) 1201,232,204 
204  IF  C Y ( J)-XiKN+l) 1206,232,201 

201  KN=KN+2 

206  UBDY*1. 

GO  TO  207 

123  L80Y3 1  • 

207  CALL  REGION 
ASSIGN  75  TO  Ll 
GO  TO  2202 

75  IF  ( IT. GT. 18000)  GO  TO  320 
755  I  =  I  + 1  NT { Y { NKYV ) )+2 

76  J3J+INT(Y(NKXV) )+2 

IF  tV.NE.O.)  GO  TO  540 

I  IST+l 

NFuNCT -NCFCT 

DO  310  NY*1 ,NYR 

NHALF=NXC/2 

DO  300  MC-ltNHALF 

IF  (Y(NFUNCT).NE.O. )G0  TO  300 

Y(NFUNCT)*YU*3) 

300  NFUNCT»NFUNCT*NYR 
NHALF-NHALF+l 
DO  305  MC3NHALF,NXC 
IF  { Y ( NFUNCT ) .NE . 0. ) GO  TO  305 
IPNPl3H>INT(Y{ I  +  l)  l  +  l 
Y ( NFUNCT ) *Y( J  PNPi 1 
305  NFUNCT  =*NFUNCT+MYR 
NFUNCT  »NCFCT*NY 
310  I  =  IMNT(Y(  !♦!  )  H-2 
540  RETURN 

320  WRITE  (6,625)  Y ( J I 
STOP 
END 
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HBFTC  WAGRDS 

THIS  ROUTINE  GETS  PSI  VALUES  OF  THE  NEAREST  POINTS  IN  FOUR 
DIRECTIONS  FROM  THE  POINT  IN  QUESTION  AND  THE  DISTANCES  TO  THESE 
POINTS 

DIMENSION  B( 1),Y{ 18000) ,SLINES(30) 

COMMON  B,NXC,NYR,NKXVCf NKYVC , NDXC ,NDYC t NCFC T , NBE T AC , K J , B I TSt 1ST ,  0 

LHONE,HTWO,ON£K,THOK, I, IN, J, JN, IKJ,XIJt JIK, IJK* I T , NDX , NDY , NFUNCT , Y,  0 

2SLINEStHl,H2,TKl, TK2,NBl,N82,NB3,NB4#NCl,NC2,NC3,NC4*Nl,N2,U0,Ul,  0 

3U3  » U2  » U4, 1 PRINr  NRSDUE  »NSL » NSLNC  »NGRDC  » NGRDXC  »NGRDYC » AXYM»NGRAD»  0 

4NGRADX»  NGRADYf  LN>CSCHf NG»  EXTRAP , SSRCH»NSWPS  t OMEGA, CVGS» CONST »  0 

5STREAM,KNTR, A,C,D,E,FfG»HKl,HKl2|HI2tHK2#BY,8X,NKXV,NKYVtN3»  0 

6P#Q»R»S»TfVfW»XC»YC» IH»NTAP6  0 

NI  =  NFUNCT+NYR  0 

N3=NFUNCT-NYR  0 

HONE=Y ( NDX )  0 

HTWO=Y { NOX— I )  0 

ONEK=Y I NDY )  0 

TWOK=Y ( NDY- 1 )  0 

UO=Y { NFUNCT )  0 

U I=Y ( Nl )  0 

U2=Y ( NFUNCT+ 1 )  0 

U3=Y { N3 )  0 

U4=Y ( NFUNCT—l )  0 

IJK*0  0 

I K J=0  0 

K I  J=0  0 

J IK=0  0 

IF  (Yl IN)-Y{ J )+Y(NOy-l) )5l,5l» 150  0 

150  HTWO  =  Y(J)-YIIN)  0 

U3=Y (  IN  +  l )  0 

I  JK  =  1  0 

51  I F { Y { IN+2 )-Y(J)-Y(NOX) ) 151,52,52  0 

I  5  I  HONE  =  Y ( IN  +  2J-YI  J)  0 

Ul-YUN+3)  0 

I K J  =  1  0 

52  IF(Y( JN)-Y( I)+Y{N0Y-1))53,53» 152  0 

152  TWOK  =  Y ( I )->YUN)  0 

U4  =  Y ( JN+  l )  0 

K  I J  =  l  0 

53  IF(Y( JN  +  2J-YI I J-Y(NOY)) L  53  r 540,540  0 

153  ONEK  =.  Y  (  JN  +  2 )-Y {  I )  0 

U2  =  Y ( JN  +  3 )  0 

J  I K=  1  0 

540  RETURN  0 

END  0 
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SUBROUTINE  ACDEFG  0 

DIMENSION  9(1) »Y( 18000) , SL INESC 30) 

COMMON  b,nxc,nyr,nkxvc»nkyvc»ndxc,ndyc»ncfct »NBETAC»KJ»BITS»l$T, 
IH0NE»HTW0,0N£K,  TWOK, I , IN , J , JN , IK J ,K I J , J IK, I JK , I T ,NDX,NOY, NFUNCT, Y, 
2SLINES,Hl,H2,TKi,TK2,NBl,N82,NB3,NB4,NCl,NC2,NC3,NC4,Nl,N2»U0,Ul» 
3U3,U2,U4, IPR IN»NRSOUE#NSL*  NSLNC  f NGRDC»NGRDXC»NGRDYC»  AXYM»NGRAO» 
4NGRADX»NGRA0Y,LN»CSCH»NG» EXTRAP, SSRCH» NS WPS » OMEGA, CVGS* CONST • 

5S  TREAM, KNTR, A,C,D»E»F,G»HK 1 , HK 12 » H12, HK2 *  B Y» BX » NKXV, NKYV* N3 * 
6NLRCAL,Q,R,S,T,V, W,XC,YC,IH,NTAP£ 

IF  ( AXYM}  1,2,1 

1  IF  |Y(I))  3,4,3 

4  C*2. 

2  E*0. 

GOTO  7 

3  E=STREAM/Y ( I ) 

7  RETURN 

ENO 


i 
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$ I8FTC  WACNST 

»  CNSTNT  SUBROUTINE  FOR  CALCULATING  CONSTANTS  OF  DIFFFERENCE  EQUATIONS 
SUBROUTINE  CNSTNT 

THIS  ROUTINE  CALCULATES  THE  WEIGHTS  TO  8E  GIVEN  THE  PSI  VALUES  OF 
THE  FOUR  SURROUNDING  POINTS  IN  THE  RELAXATION.  THESE  CONSTANTS  ARE 
STORED  BY  SEARCH  IN  THE  ARRAY  THAT  BEGINS  WITH  YINBETAC) • 

DIMENSION  B{ I),Y118000),SLINES(30) 

COMMON  B, NXC  »NYR»  NKXVC, NKYVC , NDXC » NDYC »NCFCT ,N8ETAC»KJ»8ITS,IST, 
lHONE,HTWO,ONEK,TWOK, I, IN,J, JN , IX J , K I J , JI K , I JK, IT»NDX,NDY,NFUNCT, Y, 
2SLINES,H1,H2,TK1,TK2,NBI,N62,NB3,NB4,NC1,NC2,NC3,NC4,N1,N2,U0,U1, 
3U3»U2»U4, IPRIN#NRSDUE»NSLiNSLNC»NGRDC  »NGRDXC»NGRDYCf AXYM»NGRAD, 
4NGRADX»NGRADY,LN,CSCH,NG, I XTRAP,SSRCH, NSWPS, OMEGA, CVGS, CONST, 
5STREAM,KNTR, A, C, D, E , F , G, HK 1 , HK12, H12,HK2, BY , BX, NKXV,NKYV,N3, 
6P,Q,R,S,T,V,W,XC,YC, IH,  NTAP6 
COMMON  U8DY,LB0Y 

THE  CONSTANTS  I KJ , J I K, I JK, K I J  ARE  SET  EQUAL  TO  l  IF  LESS  THAN  A 
MESH  DISTANCE  AWAY  FROM  A  80UNDARY  POINT. IN  THIS  CASE  THE  FUNCTION 
VALUE  AT  THE 

BOUNDARY  IS  INCLUDED  IN  THE  CONSTANT  TERM  YIIT+4),  AND  THE 
CORRESPONDING  CONSTANT  TERM  IS  SET  EQUAL  TO  ZERO. 

222  BA= ( HTWQ»D+2. » A ) / ( HONE* ( HONE+HTWOI ) 

8B=  (  2.  »C  +  TWOK*£  )  /  { ONEKU ONEK+TWOK ) ) 

BC*  (  2.*A-H0NE*0 )  /( HTWO*  (  HDNE+HTWO ) ) 

80=*  (  2 .*C-ON£K*E  )  /  {  TWOK*  I ONEK+TWQK ) ) 

BE=-F*BA+BB+BC+BO 
8E=»  l./BE 
Y ( I T+4 ) *-G*8E 
I F  C  IKJ)202,20l,202 

202  BA  *  BA*Y ( IN+S ) *8E 
Y(IT)*0. 

Y(  IT+4)a8A+Y{ I T+4 1 
GO  TO  210 
201  Y { I T ) =8A*BE 

210  I  F { JIK)204,203,204 

204  B3  =  BB*Y{ JN+3) *3E 
Y  ( I TU  )  =0. 

Y { I T+4 ) SBB+Y ( IT*  4) 

GO  TO  211 

203  Y (  ITU)  *BB*BE 

211  IFUJK)206,205,206 

206  8C=BC*YUNU)*BE 
Y (  I T+2 ) =0 • 

Y( IT*4)=BC+Y( 1T*4) 

GO  TO  212- 

205  Y (  IT*2)*9C*B£ 

212  IFC  ,<IJ  )  208, 20T,  208 
208  BD=30*Y(JNU)*B£ 

Y(  I T+3 ) *0. 

*  YUT  +  4)*BD*YUT*4) 

GO  TO  213 

207  Y { I T+3 ) *BD*8£ 

C  ANY  DERIVATIVE  ENFORCEMENT  BEGINS  WITH  STATEMENT  213 
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213  IF  (IXTRAP)  217,227,217 
217  IF  (UBOYI  1,2,1 

1  Y{ IT+3)=Y(IT+3)*Y{ ITU) 
Y(ITU)*0. 

GO  TO  225 

2  IF  (LBDYJ  3,227,3 

3  Y(  I T+ l } =Y {  IT  +  3M-Y(  IT*1) 
Y ( IT  +  3)*0» 

225  U80Y=*0. 

LBOY-O 
227  IT*lT+5 
RETURN 
END 


%> 
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SI8FTC  HASHEP 

SUBROUTINE  SWEEP? PUNT,. NSWEEP, HEDGE, FllH) 

DIMENSION  STORE! 25) 

DIMENSION  B ( I ) , V 1 18000) , SL INES t  30) 

COMMON  B,NXC,NYR,NXXVC,NKYVC,NDXC,NDYC,NCFCT,NBETAC,KJ,BITS,IST, 
1H0NE , HTWO, ONEK, TWOK, I»!N,J,JN,IKJ,K/J»JIK,IJK,I T, NDX , NDY, NFUNCT, Y, 
2SLlNES,Hi,H2,TKi,TK2,NBl,NB2,NB3,NB4,NCl,NC2,NC3,NC4,Nl,N2,.U0,Ui* 
3U3»U2»U4, IPRIN, NRSDUE  »NSL»NSLNC, NGRDC, NGR0XC»NGRDYC, AXYMf NGRAD, 
4NGRA0X »NGRAOY,LN,CSCH,NG, EX TRAP, SSRCH,NSWPS, OMEGA, CVGS, CONST, 
5STREAM, KNTR »A,C»D,E,F,G»HK1,HK12»H12,HK2,BY»  8X»NKXV»  NXYV, N3 , 
6P,Q,R,S,T,V,W,XC,YC,IH,NTAPE 
COMMON  /VEL/  NGRDVC 

253  FORMAT  ( / 13H  SHEEP  NUM8ERI5,  9H  RES  I DUE*E15, 7, 

I 10H  AT  COLUMN, 14, 4H  R0H,I4,25H  WITH  CONVERGENCE  HIST0RY//I6F18.6) J 
255  FORMAT  (36H0PR0GRAM  LEFT  SUBROUTINE  SWEEP  AFTER, I8,12H  SWEEPS  WITH 
XF12.6»8H  RESIDUE//) 

686  FORMAT  (IX, 4110) 

750  FORMAT  ?1H1) 

WRITE  16,750) 

IN*0 

REWIND  8 
NSKEP1*NSWEEP 

WRITE! 8) ( Y  C INM), INM*NBETAC, 18000) 

316  MK=NXC»NYR 
Nsl 

310  DO  311  NNN* 1,25 

311  STORE (NNN) *0# 

MN*2 

DO  112  JNl=*i,NSW£Pl 
IN* IN+l 
JN*JNl/2+l 
MN=3-MN 

204  GOTO! 1,2 ) ,MN 

1  NBETA=NBETAC 
NADD* l 

NFUNCT*NCFCT 
GO  TO  3 

2  NBET A=N0E T A- 1 
NADD=-i 

NFUNCT=NFUNCT~l 

3  P.SD1*0, 

DO  12  1 1 3 1 ,NXC 

DO  12  1=1, NYR 

IF?Y(N8ETA)-BITS)4,10,4 

4  TEMP=YINFUNCT) 

GO  TO  ( 6, 5 ) , MN 

5  NBETA=NBETA-4 

6  TEMP1=RELAX(NBETA, NFUNCT) 

I F (  (I.EO.l).OR.n.GE.NYR))  GO  TO  601 
RSD=TEMP1-TEMP 
19  RSD=R$D*OMEGA 
17  Y(NFUNCn=TEMP+RSO 
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601  GO  TO  ( 7, 8 ) *  MN 

7  NBETA=NBETA+4 

8  IF! JN-1)9»9, 10 

9  IF  (RSD1-ABS  I RSO ) )  99,10,10 
99  RSD1=*A8S  (RSO) 

R$D3=RSD 
RS04-RS03 
NLOC=NFUNCT 
NL0C1=N6ETA 
GOTO( 20,21) »MN 

20  NL0Ci=NL0Cl-4 

21  NFU=(NFUNCT-NCFCT+l)/NYR 
NFUN=(NFUNCT~NCFCTU)-NFU*NYR 
IFINFUNJ97, 98,97 

98  NFUN=NYR 
97  NFU=NFU+l 

10  NBETA-NBETA+NAOO 
12  NFUNCT=NFUNCT*NADD 
400  GOTO{ 301, 1 12 ) #MN 

301  STOR£(JN)3Y{ NLOC ) 

GO  TO ( 1 12 , 299 ) ,N 

299  IF(JN“2)302, 112,303 

302  IF(RSD1~CVGS)502»304*304 

303  I F ( (STORE  (  JN)-STORE ( JN-l } }/(  STORE  ( JN-U -STORE  { JN-2  >? 1298,112, U2 
298  N'2 

GOTO  306 

304  IENO=NCFCT*MK 

WRITE  (8)  (Y(IJN), IJN=NCFCT, IENO) 

318  BACKSPACE  8 
112  CONTINUE 

GO  T0( 306, 320) »N 

320  F AC T= ( STORE ( 6) -STORE ( 5) ) /( STORE ( 2) -STORE ( l ) ) 

I F ( FACT- . 98)329, 330, 330 

329  FACT=( 1 . / ( l*-FACT)-l»  )* HEDGE *!• 

F  ACT  =  AMI N 1 ( F ACT , FL IM  ) 

GOTO  331 

330  FACT=FLIM 

331  NBETAX*NBET  A+i-NCFCT 
I ENO=NBETAC+MK 

READ  (0)  (Y( IJN) , IJN-NBETAC, IENO) 

322  NFUNCT=NCFCT 
NBETA-NBETAC 
DO  321  I » 1 , MK 

Y(NFUNCT)*(Y(  NFUNCT  ) -Y(  NBETA ))  »FACT  *-Y  ( N8ETA ) 

NFUNCT3NFUNCTM 
321  N8ET  A=NB6TA+1 
324  REWIND  8 

IENO=NBETAC+MK 

REAO  (8)  (Y( IJN), IJN-NBETAC, IENO) 

306  WRITE  (6,253)  IN,RSD3,  NFU, NFUN , t STORE ( I ) , !*i , JN) 

16  IF(  IN-NSWPS.'  307,  502,502 

307  GQTO(308,309),N 
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303 


309 


502 

105 

61 


NSWEPl=NSHEEP 
N*  2 

GOTO  310 

NSWEPl*ll 

N-i 

GOTO  310 

WRITE  (6,255)  IN,RSQ4 
00  105  i=NGROXC,NGROY€ 
YU 

RETURN 

END 
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000000 


tIBFTC  WARELX 

FUNCT IGN  RELAX! N8»NFUNCT)  0 

DIMENSION  B ( 1 ) » Y 1 18000) »SLINES(30) 

COMMON  B, NXC  *NYR ,NKXVC , NKYVC  , NDXC-*N0YC  »NCFCT»N8ETAC  »K J  »BITS»IST» 
1H0NE*HTW0,0NEK,TW0K, t,IN,J,JN, iK J,K IJ , J I K, I JK, I T , NDX»NDY,NFUNCT, Y* 
2SLlNES,Hl,H2,TKltTK2,NBl,N82,NB3,NB4,NCl,NC2,NC3,NC4,Nl,N2,U0,Ul, 
3U3,U2,U4, IPRINf  NRSDUEf NSL » NSLNC»NGRDC  *NGRDXC»NGRDYC* AXYM»  NORADf 
4NGR ADX,NGRADY,LN,CSCH,NG, EXTRAP, SSRCH,NSWP$, OMEGA, CVGS, CONST, 
5STR£AM,KNTR, A,C»D,E,F,G»HKL,HKl2,Hi2,HK2»BY»BX»NXXV*NKYV, N3» 
6P,Q,R,S,T,V,W,XC,YC, IH,NTAPE 
N1  =  NFUNCT  +NYR 
N3=NFUNCT-NYR 

^  RELAX=Y(Nn*Y(NBH-Y{NFUNCm)*Y{NBUH’YlN3l»Y{NB*2)«‘ 
lY(NFUNCT-l)*Y(N8+3)<-YINB4-4) 

RETURN 
END 
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% LBr TC  WALNK2 

SUBROUTINE  VELOC 

DIMENSION  8(  1)  ,-Y  (18000).,  SL INESI  3(71' 

COMMON  B,NXC,.NYR,.NKXVC*NKYVC,NDXC».NDYC«NCFCT»N8eTAC,KJ,8ITS,IST, 
lHONE,HTWO,ON£K,TWOK,  1,-IN,  J:,  JN,IKJ,KIJ*  JIX,I  JK,.I  T,NDX,NDY,NFUNCT,Y,. 
2SLINES,HI,H2frTKUTK2,.N8ltN82»N83>NB4fNCl,NC2»NC3»NC4,Nl,N2.UO.Ul* 
3U3,U2,U4,  IPRlNr.NRSDUE,N$L,NSLNC,NGRDC»NGRDXC*NGRDYC»AXYM,NGRAO,. 
4NGRADX  ,NGRADY#LN».CSCKfNG»ND£R  IV,  SSRCH,NS'WPS,  OMEGA  ».CVG5»CQN$r  r 
5S.TR£AM,KNTR»  A.»C »  0>E»F».G»HK<l,HK12».H22,HK2,8Y ,  BX , NKXV»NKYV»N3» 
6P»Q»R»S,T»V»W»XC  # YC » IH»  NTAPE 

COMMON  /COM2/  X$CL,Y9CL, PUNCH,  I  TIE , YO, Yl , XL l , XL2 » YL 1 , YL2 • XSCLA, 
1YSCLA,XL1A,YL1A,XL2A,YL2A 
COMMON  /VEL/  NGROVC 
DIMENSION  ITL6I24) 

650  FORMAT  { 1H1, 3X,2HX^F14.7//I 

651  FORMAT  { 11X, 1HY, 16X , 3HPS I ,12X , 12HX-DER I VAT  I V£,.6X, 12HY-DERI VATI VE# 
19X , 8HVELCC ITY// ) 

652  FORMAT  ILX,.5F18.7) 

686  FORMAT  (IX, 4110} 

ZER0=O. 

7  V*t. 

N2=»l 

112  DO  1  NN*1, 2 
NGRACeNGROC 
CALL  ICMPGS 
1  CSCH=C$CHU. 

16  K*=NGROYC 
J  ~K  J 

NFUNCT=.NCFCT 

NGRADX=NGROXC 

NGRACY=NGROYC 

NGRACV=NGRDVC 

CO  250  KP= l , NXC 

LINE=0 

WRITE  (6,650).  Y(J) 

WRITE  {6,t>5U 
i - i sr* i 

DO  251  KPX=UNYR 

•IF  (Y(NGRADX)  «NE,8tT31  GO  TO  200 
WRITE  (6,652)  Y( I) , ZERO, ZERO, ZERO»ZERO 
LCNE*UNEf  1 

IF  (LINE. GT„ 52)  GO  TO  101 
GO  TO  201 

200  WRITE  (6,652)  Y( I) ,Y (NFUNC  T),Y(NGRADX),Y(NGRADY) ,Y(NGRAOVS 
LIN6=LIN£+1 

IF  ( L INE.GT , 52 1  GO  TO  101 

201  NFUNC  T  =.  NFUNC  T*1 

T  s  I  + 1 N  T  ( Y  ( I ♦ l } )  ♦  2 
NGRACX=iNGRAOX#-i 
NGRACY=NGRAQY*i 
251  NGRACV^NGRAOVM 
250  J*J+INT(Y(  J+U  )*2 
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331  RETURN 
iOi  L  CNE=0 

WRITE  (6,650)  ¥U) 
WRITE  (6,651J 
GO  TO  201 
END 
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01 


H8FTC  WAICMP 

SUBROUTINE  ICHPGS 

C-  THIS  ROUTINE  IS  CALLED  TWICE*  THE  FIRST  TIME  IT  CALLS  GROST1  AND 
C  D6RIV  FOR  INTERIOR  POINTS.  THE  SECOND  TIME  IT  CALLS  SUBBi  AND 

C  OERIV  FOR  BOUNDARY  °OINTS. 

DIMENSION  BU),Y{  10OOO)  ,SLINESI30i 

COMMON  B»NXC»NYR»NKXVCf NKYVC#NDXC*NDYC  »NCFCT »NBETAC»KJ»8ITS»I$T<> 
1H0NE,HTW0,0NEK,TW0K, I , IN»J.JN,IXJ#KI J, JIK,IJK,IT»NDX,NDYtNFUNCT|Y, 
2SLINES,Hl,H2,TKl,TK2*NBi,NB2,NB3,N84,NCi,NC2,NC3,NC4,Nl,N2»U0,Ul, 
3U3»U2»U4, IPRIN,NRSDUE, NSL, NSLNC»NGRDC, NGRDXC»NGRDYC, AXYM,NGRAD, 
4NGRADX, NGP  iHY, LN,CSCH,NG, EXTRAP, SSRCH,NSWPS, OMEGA vCVGS, CONST, 
5$TR£AM,KNTR,  A,C,D,E,FrG,HKWHXI2fH12,HK2#6Y,6X,NKXV,NKYY»N3» 
6P»Q,R»S,TfV,W»XC,YC»IH*JH,KH 
COMMON  /VEL/  NGRDVC 

*  INITIALIZE  ALL  SUBSCRIPTS  TO  PICK  UP  INFORMATION  IN  Y  ARRAY. 

• 

NCI*  l 

NGRADX*NGRDXC-l 

NGRAOY=NGRDYC~l 

NGRADV*NGRDVC-I 

NKXV’*NKXVC-1 

NFUNCT*NCFCT-i 

NDX-NDXC”! 

J*K  J 

» 

#  SET  LOOP  TO  PROGRESS  FROM  COLUMN  TO  COLUMN. 

* 

DO  76  NX*1 » NXC 

*  CHANGE  SUBSCRIPTS  AFFECTED  BY  A  CHANGE  OF  COLUMN. 

» 

NDY=NDYC-l 
NKYVXNKYVC- l 
NDX=N0X+1 
I  =  IST+I 
NKXV--NKXV+I 

»  SET  LGOP  fO  PROGRESS  ALONG  A  COLUMN. 

DO  75  NY=l,NYR 

*  CHANGE  SUBSCRIPTS  AFFECTED  BY  A  CHANGE  ALONG  A  COLUMN 

NGRAOX-NGRADX*! 

NGRADY=NGRADY*l 
NFUNCT=NFUNCT+l 
NGRADV*NGRA0V+1 
NKYV  =  NKYV«-l 
NOY=NDY+ l 

*  SET  SUBSCRIPTS  TO  FIND  BOUNDARY  INFORMATION  FROM  THE  INPUT  ARRAY. 

* 

JN*J>2 

IN*I+2 

KYVA*Y { NXYV } /4. 

KYVB  *  KYVA 
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*  LOOP  10  ANO  112  AR£  NECESSARY  FG*  COMPLEX  BOUNDARIES* 

DO  10  IS31»KYV.A 

e 

*  LOCATE  POSITION  OF  POINT  IN  THE  ENTIRE  SYSTEM* 

* 

•  THE  POINT  HILL  EITHER  BE  INTERIOR, EXTERIOR  0R60UNDARY* 

» 

IF{Y(IN)—YCJJ)l,2»3 

1  I F ( V C IN+2J-YIJ)  }  110,22,7 
7  KXVA=Y(NKXV)/4. 

KXV8=KXVA 

DO  112  J$= 1 »KXVA 

I F ( Y{  JNI-Y { I JU1, 222, 3 

11  I F { Y C  JN+21-YU)  )12V,2222,13 

• 

*  THE  DERIVATIVES  FOR  INTERIOR  POINTS  ARE  COMPUTED  ON  THE  FIRST  SWEEP. 

13  IF  (CSCH-1.)  113,75,75 
113  CALL  GRDST1 ( NXA,NYA) 

GO  TO  90 

12  KXVB3KXVB-1 

IF  ( KXVB )  3,3,112 
.  112  JN3 JN+4 

110  KYVB  «  KYVB-1 

IF  (KYVB)  3,3,10 
10  IN=*  IN+4 
GO  TO  75 

♦ 

*  THE  DERIVATIVES  ON  BOUNDARIES  ARE  COMPUTED  ON  THE  SECOND  SWEEP. 

• 

2  N3l3i 
GOTO  202 

22  NB 1=2 
GO  TO  202 
222  NB1=3 

GO  TO  202 
2222  N81---4 

20 2  IF  (CSCH-1. 1  75,31,75  ' 

# 

*  SUBROUTINE  SU881  SUPPLIES  THE  INFORMATION  FROM 

* 

*  WHICH  THE  BOUNDARY  DERIVATIVES  ARE  COMPUTED. 

31  CALL  SUB81 ( NXA»NYA) 

* 

•DERIVATIVES  ARE  COMPUTED  USING  CENTRAL, FORWARD  OR  BACKWARD  DIFFERENCING 

♦ 

90  CALL  DERIV  { NXA , NYA , NXYAA, BX , 8Y , G3 ) 

» 

*  ON  RETURN  FROM  DERIV, BX  AND  BY  ARE  THE  X  AND  Y  DERIVATIVES  . 

Y ( NGRADX } *BX 
Y (NGRADY )*BY 
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Y  J  H^RADV ) *G3 
GO  TO 

3  Y {NGRAOX)*BITf 
Y { NGRADY ) *81TS 

BITS  ARE  STORED  IN  THE  X  ANO  Y  DERIVATIVE  ARRAY  FOR  EXTRAPOLATION  USES 

75  I3I+INT(Y{ NKYV ) ) +2 

76  J=J+I NT { Y ( NKXV ) )*2 
2.6  RETURN 

ENO 


1073 
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$  I BFTC  WAGROi 

SUBROUTINE  GRDSTKNXA.NYAl 
DIMENSION  B(i),Y{  18000)  »SUNESC30) 

COMMON  8,NXC,NYR,NKXVC,NKYVC,NDXC,NDYC,NCFCT,NB£TAC,KJ,BITS,ISTt 
1H0NE,HTW0,0NEK.TW0K» I,IN,J, JN, IKJ ,K l J , JI K, I JK, I T, NDX, NDY.NFUNCT, Y, 
2SlINES,Hl,H2,TKi,TK2,N81,N82,NB3,N84,NCl,NC2,NC3,NC4,Nl»N2,U0,Ul, 
3U3.U2.U4, IPRIN.NRSDUE.NSL.NSLNC.NGRDC.NGROXC.NGRDYC, AXYH.NGRAD. 
4NGRA0X, NGRAOY.LN.CSCH.NGfEXTRAP.SSRCH, NSWPS, OMEGA, CVGS. CONST, 
5STREAH,KNTR,A,C,D,E,F,G.HK1IHKI2,H12,HK2,8Y,BX.NKXV,NKYV,N3, 
6P»0,R,S,T»V.W,XC»YCfIH 
NXA-l 
NYA*l 

NI  *  NFUNCT+NYR  . 

N3aNFUNCT-NYR 
Hl*Y  (NDX) 

H2*Y ( NOX-1 I 
TK1=Y (NDY) 

TK2=Y { NDY-I ) 

UOaY I NFUNCT ) 

UI-Y(Nl) 

U2=Y(NFUNCT*U 

U3-YIN3) 

U4*Y  t  NFUNCT- 1 ) 

I  JKaO 
IKJ*0 
K!  J30 
J  IK30 

IF  (Y(IN)-Yt  J)+Y(NOX-U  i 51,51  ei50 

150  H2  3  YIJ)-Y(IN) 

U3-Y ( IN+ I ) 

I  JK=l 

51  I F  ( Y  (  IN  +  2)”*Y{J)-Y(N0X) )  151,52,52 

151  HI  =  Y ( IN+2J-YI J) 

Ul-Yl IN+3) 

I  K  J  *  I 

52  I F { Y C JN ) — Y { I )+YINDY-i)  153,53,  152 

152  TK2  =s  Y(I)-Y(JN) 

U4-Y( JN+l) 

K  I  J3 1 

53  IF{YtJN+2)-Y{ I)-Y(NOY)) 153,540,540 

153  TKi  *  Y{JN+2)-Y(I) 

U23Y( JN+31 

JIK=l 

.  5*0  RETURN 
END 


$1 BFTC  WASUBi 
*SU8B1  00 

SUBROUTINE  SUB3UNXA,NYA) 

DIMENSION  6il),Y« 13000) ,SLINES( 30 ) 

•  COMMON  b,nxc,nyr,nkxvc,nkyvc,ndxc»ndyc»ncfct»nbetac,«j,bit$»ist. 

IHONE ,  HT WO,  ONEK ,  T WOK  ,  I ,  I N » J  ,  JN ,  I  < J  ,K I J  ,  J I K ,  I  JK ,  I  T » NDX , NO* , NFUNCT »  Y , 
2SLlNES,Hl,H2,TKl,TX2,NBl,NB2,N83,NB4*NCI,NC2,NC3,NC4fNlsN2JU0,Ul, 
3U3,U2,U4, IPRIN,NRSDUE,NSL,NSLNC,NGRDC,NGRDXC,NGRDYC»AXYM,NGRAD, 
4NGRADX,NGRADY»LN,CSCH,NGi EX TRAP, $ SRCH, NS WPS »GKEGA»CVG$, CONST  > 
5STREAM,KNTR,A,C,D,£,F,G,HKl,HKi2,Hl2,HK2,BY,BX?NKXV,NKYV,N3t 
6P,Q,R,S,T,V,W,XC.YC»IH 
301  FORMAT  (28H  THERE  IS  GARBAGE  IN  MEMORY//) 

*  THIS  SUBROUTINE  FINDS  DATA  TO  COMPUTE  DERIVATIVES  FOR  BOUNDARIES. 

*  CHOOSE  WHICH  BOUNDARY  IS  BEING  WORKED  ON. 

MMH=*3 

N1*NFUNCT+NYR 

GO  TO  1 12, 13, 14, 15) , NBl 

* 

*  LEFT  BOUNDARIES  ARE  IN  QUESTION. 

* 

12  NP*NGRAOY+NYR 
N3*N1+NYR 
UO*YUNM) 

U1=Y ( N 1 ) 

Hi*YINDX) 

NXA*2 

IF  (YtN3)“BITS)30,2»l 

1  U3=Y ( N3 I  v 

H2*Y(NDX*i) 

GO  TG  100 

2  U3=Y(IN»3) 

H2  =  ABS  (YI  JJ-YUN+ZH-Hl 
GO  TO  100 

*  RIGHT  BOUNDARIES  ARE  IN  QUESTION. 

• 

13  NP*NGRADY-NYR 
NXA*3 

N 1=NFUNCT~NYR 
N3=N1-NYR 
UOaY { IN+3 ) 

U i=Y { Nl ) 

Hl=Y{NDX-i) 

IF  (Y(N3)~81T$)30, 20,31 
3 1  U3-YIN3) 

H23Y{ NDX-2) 

GO  TO  100 
20  U3  =  YUN+1) 

H2  =  ABS  ( Y C J)-Y( IN) )-Hl 

« 

*  SET  UP  DATA  TO  FACILITATE  FINDING  Y 'DERIVATIVES 
100  T5=»Y(N0Y~2) 
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T6*YtNDYUJ 
U2*Y(NFUNCT*il 
U6  =  Y ( NFUNCT  *21 
U4=Y(NFUNCr»i: 

U5  =  Y I NFUNCr-2  2 
TKl=Y{NDY ) 

TK2sYCN0Y-l ) 

* 

*  IS  THE  BOUNDARY  POINT  A  CORNER  OF 


THE  CIRCUMSCRIBED  RECTANGULAR  REGION 


I F( NOYC-NDY )  22*10*22 

22  I F(NDYC~NDYfNYR- 1)23*7# 23 

23  IF(U2-B1TS)30*7»6 

6  IF(U4-BITS)30tl0,90 

* 

»  THE  BOUNOARY  POINT  IS  A  TOP  CORNER. 


7  IF(U4-SITS)30,8,U 
11  IF(U5-BIT$)30,8*2t 
21  U2-U4 
U4*U5 
TK1*TK2 
TK2=*T5 
NYA*3 
GO  TO  9 

» 

*  THE  POINT  I  SA  BOTTOM  CORNER  OF  THE  RECTANGLE 

« 

10  IF  { U2-B I TS ) 30» 8 , 18 
18  IF  (U6-BITS)3Q,8, 120 
120  U4=U6 
TK2=T6 
NYA=2 
GO  TO  9 

# 

.»  A  LINEAR  EXTRAPOLATION  IS  MADE  WHEN  DATA  IS  INSUFFICIENT. 


8  GO  TO  <112*113* 114* 115) *NBl 

112  NP1=NP+NYR 
D1  =  H1 
D2=H1*H2 
GO  TO  116 

113  NP1=NP-NYR 
Dl  =  H2 
D2*HUH2 

116  NYAa4 

BY  =  YINP1)-(D2/Di  )*  ( YINP  U-YINPH 
GO  TO  9 

114  NPl=NP*l 
Ol-TKl 
02*TKUTK2 
GO  TO  117 
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115  NP1*NP-1 
D1*TK1 
D2=TKi+TK2 
117  NXA*4 

8X*Y ( NP l ) -( D2/D1 )*  S Y(NP 1 ) -Y {NP ) ) 
GO  TO  9 

14  N3=NFUNCT-NYR 
NYA  *2 
U0*Y(JN+1) 

U2*YCNFUNCT  +  U 
TK1*Y { NOY ) 

IF  C  Y {  NFUNCT+2J-BITS )30t 82»81 

81  U4  *  Y(NFUNCT+21 
TK2^YCNDY+1> 

GO  TO  1000 

82  U4*YIJN+3> 

TK23ABS  lYCD-Yt JN+2D-TKI 
1000  IF(Y(N3)-BIT$)30,8t30Q 
300  U3*Y { N3 ) 

H2=Y { NOX-1 ) 

IF  (YtNl)-BITS)30,8t70 
70  U1«Y(N1I 
Hi*Y (NOX) 

GO  TO  91 

15  N3=*NFUNCT-NYR 
NYA*3 

U0=Y( JN+3) 

U2*Y( NFUNCT-i ) 

TKl*Y(NDY-l) 

I F { Y { NFUNCT-2 )-B I TS ) 30»  42t  41 

41  U4  =  Y { NFUNCT-2 ) 

TK2=Y ( NOY-2 ) 

GO  TO  110 

42  U4=Y(JN+1) 

TX2  =  ABS  ( Y(  I  >-Y(  JN)  )**TK1 
110  IF  (Y(N3)-BITS)30,8,43 

43  U3*Y(N3) 

H2=Y(NDX-n 

IF  CY {Ni ) — BITS3  30»8*47 
47  U i*Y { N 1) 

H1=»Y  I NDX  J 
91  NXA*l 
GO  TO  9 
90  NYA=1 
GO  TO  9 

30  WRITE  (6»30U 

GO  TO  (8l,821.MMM 
9  RETURN 
ENO 
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SIBFTC  WAOERV 

SUBROUTINE  DERI V  INXA,NYA, NXYAA, XDER IV, YDERI V,G3) 

SUBROUTINE  DERIV  CALCULATES  DERIVATIVES  0?  PSI  WITH  RESPECT  TO  X 
ANO  Y,  AND  VELOCITIES. 

COMMON  B»NXC  »NYR»NKXVC*NKYVC»NDXC»NDYC»NCFCTfNBETAC*KJ»8lTS,  TST, 
1H0N£,HTW0,0NEK,TW0K, I,IN,J,JN,IKJ,KIJ,JIK,IJK,IT,NOX»NDY,NFUNCr,Y, 
2SlINES,HlfH2,TKl,TK2,NBl,NB2,NB3,NB4,NCl,NC2,NC3,NC4,Nl,N2,U0,Ul, 
3U3,U2,U4,IPRIN,NRSDUE»NSL,NSINC,NGRDC,NGRDXC*NGRDYC,AXYM,NGRADP 
4NGRADX,NGRADY,LN,CSCH,NG,NDERIV, SSRCH, NSWPS, OMEGA, CVGS*CONST, 
5STREAM,KNTR, A,C,D,E,F,G,HK1,HK12,H12,HK2,BY,BX,NKXV,NKYV,N3» 
6NLRCAL,Q,R,$,T,V,W,XC,YC,IH,NTAPE 
DIMENSION  B{ 1 ) * Y C  18000) , SlINESI 30) 

600  FORMAT  (1X»6F14.6) 

*  IF  NXA=1  THE  X-DERIVATIVE  IS  FOUND  BY  CENTRAL  DIFFERENCING 

*  IF  NXA=2  THE  X-DERIVATIVE  IS  FOUND  BY  FORWARD  DIFFERENCING 

*  IF  NXA*3  THE  X-DERIVATIVE  IS  FOUND  BY  BACKWARD  DIFFERENCING 

*  THE  ABOVE  CONVENTION  HOLDS  FOR  NY A  AND  Y-DERIVATIVES 

*  IF  NXA=4  NO  X-DERIVATIVE  IS  FOUND. 

FCT=l. 

IF  (AXYM.EQ.l.)  FCT *Y( I ) 

NO*  1 

100  DELTA*1. 

FG31. 

GO  TO  (3,2, 1,7),NXA 

1  OELTA*-!. 

2  H2=H1*H2 
FG=-i. 

3  GX=DEITA*(H2**2*(U1-U0)+  Hl**2* ( U0-U3) ) / t Hl*H2* (H2*FG*H1 ) ) 

4  GO  TO  ( 6 « 10) » ND 

6  XDERI V*GX 

IF  (NYA-4)  7,15,7 

7  ND*2 

NXA=NYA  • 

Ul=U2 
U3  =  U4 
Hi-TKl 
H2*TK2 
GO  TO  100 

*  G3=GRADI6NT,YOERIV=Y  ER  CROSS  DERIVATIVE »XDERI V=X-DERIVATl V£ 

10  YDERI  V.^GX 

15  CONTINUE 

G33SQRT.  IIXDERIV/FCT)**2  +  IYDERIV/FCT)»*2) 

13  RETURN 
END 
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000000 


^  ilBFTC  WALNK3 
'  SUBROUTINE  LINES 

COMMON  8,NXC,.NYR,.NKXVC*NKYVC,.NDXC,NDYC,NCFCT,NBETAC,KJ,81TStIST, 
lHONE,.HTWO,ONEX,TWOK,  l,  IN,  J,  JN,IKJ,KU,  JIK„l  JK,  !  T,NDX,NDY,NFUNCT,Y, 
2SL INES , HI » H2  ^TKl ^TK2 ►N3 l,NB2»NB3,N84,NCi»NC2»NC3»NC4,Ni»N2»UO»Ul* 
3U3,U2,U4,-I PRIN^NRSDUE ►NSL,NSLNC , NGROC  »NGRDXC ,NGRDYC , AXYM,  NGRAD,. 
4NGRADX,NGRADY,LN,.CSCH,NG,NDER I V,$SRCH, NSWPS, OMEGA, CVGS, CONST, 
5STREAM , KNTR »  A,C »  D,  E » F ».G , HKl  ,HK  12 ».H12  »HK2  »BY  ,BX  »NKXV , NKYV»N3* 
6PtQ,R,S,T»V,W,XC».YC» IH,.NTAPE,  T I  TLEl ,NCHAIN, MONTH 
COMMON  /COM2/  XSCL  ».YSCL » PUNCH,  I  TIE , YO, Y1 , XL UXL2 , YL l , YL2 , XSCLA, 

1 YSCLA , XL IA, YL1A, XL2A, YL2A 
COMMON  /CCM3/  YMAX 
COMMON/ 5-RSCOF/  FlRST,.VO,CO 
COMMON  /VEL/  NGROVC,  MR 
DIMENSION  18000} ,SLINES(30J,  ITLEI24) 

DIMENSION  B(U 
501  FORMAT  (6E12.T6I 
686  FORMAT  (1X,4U0) 

2  :ILN-0 
NNLNSM 

34  IF  ( ILN.EQ.NNLNS)  RETURN 
36  READ  .(5,501)  ANSL  ►CSCH,NNLNS ,CQ  ►VO,  AIR 
IF  (CSCH.EQ.2. i  NGRDC=NGRDXC 
IF  (NNLNSiLT «!1 ).  NNLNS*l 
NSL-ANSL 

READ  (5,501)  I SL INES ( INK), INK*i »N9L ) 

,  12  NSLNCkNGRDC 

DO  61  NG= 1 , NSL 
LN=NSLNC 
11,1  N2*  l 

CALL  MARC 
1111  NGRAOfcNGRDC 

60  CALL  SORT  (  INT(CSCH)-l); 

61  CONTINUE 
27  I LN= ILN+ 1 

GO  TG  34 
ENO 


1! 

If 

If 

If 

1! 

li 

If 


if: 


A  . 

I! 

I! 

if 


If 

If 

1< 
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SIBFTC  HASGSC 
•SGSRCH 

SUBROUTINE  SGSRCH 

DIMENSION  Y( 18000) ,5LINES{30) 

COMMON  B,NXC,NYR, NKXVC,NKYVC,NDXC»NDYC,NCFCTf NBETAC, KJ»8IT$» 1ST* 
1H0NE»HTW0,0N£K,TW0K» I , IN* J , JN *  I X J *K I J * JI K, I JK, I T ,NDX,NDY,NFUNCT, Y* 
2SLINES,Hl,H2,TKl,TK2,NBl,N82,N83,N84,NCl,NC2,NC3,NC4,Nl,N2,U0,UW 
3U3»U2»U4» IPRIN,NRSDUE,NSL,NSLNC,NGRDC,NGRDXC»NGRDYC,AXYM,NGRAD, 
4NGRADX,NGRADY»LN,CSCH,NG»  EX TRAP*  S SRCH*NSHPS» OMEGA *CVGS»CONST * 
5STREAM,KNTR»A,C»D»E»F,G»HK1,HKI2»H12»HK2  *  BY*BX»NKXV»NKYV*N3  * 
6P»Q.R*S*T»V»W*XC»YC» IH,NTAPE 
COMMON  /ZCOM/  Z 
COMMON  /LNK5/  NV 
COMMON  /COM/  U7* Hi 
C0MM0N/3DRY/  JB 
COMMQN/PRSCOF/  FIRST»VO 

600  FORMAT  QX,F14.7) 

601  FORMAT  I48H  STORAGE  OVERRUN  WHILE  CALCULATING  STREAM  LINES  ) 

602  FORMAT  (1X,5F14.5) 

N=NG 

IF  ( CSCH-3« )  57,81,82 

81  Ni=NGRADX+NYR 
N3=NGRA0X-NYR 
UO=Y ( NGRADX I 
Ul*Y(NU 
U2*Y (NGRADX*!) 

U3XY ( N3 ) 

U4=Y { NGRADX-l ) 

\  GO  TO  57 

'  82  Nl=NGRADY+NYR 

N3=NGRADY-NYR 
UO=Y ( NGRADY ) 

U IaY { N1 ) 

U2  =  Y  ( NGRADY+1 ) 

U3-Y(N3) 

U4- Y ( NGRADY- 1 ) 

57  N2-LN 
V=0. 

IF  { K I J )  84,85,84 

84  IF  (CSCH-3.)  85,2,2 

85  U6-U4 
U73U2 

W I*0N6K 
W=-TWOK 

CALL  SUB58  (U6,N> 

IF  (VI  3,2,3 
3  Y { LN  + 1 ) XY  t  J ) 

Y ( LN+2 ) -Y (I ) +XC 

IF  (KIJ.EQ.1.AND.J8.EQ.1)  GO  TO  31 
NVL=NVU 

Y  { LN*3  ) 3Y { NV ) •*■(  Y { NV i  ) -  Y( NV ) ) »XC/W 
GO  TO  32 
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31  W*TWQK 

YaN+3)*{W+W!)*lYtNV)-Y(NV-»-in/Wl+YCNV+n 

32  Y(LN+4)*i*-{Y(LN+3)/VQ)**2 
LN*LN+5 

2  V*0o  l 

12  Y*0*  i 

IF  (CSCH-3.)  69*90,90 
89  IF  C  J  IK  I  22 *T3»22 

22  H*0NEK  ‘ 

wi^-rwox 

Ufe«U2  i 

U7*U4 

CALL  3U858  (U6tN) 

IF  (V)  21,73,21  '  * 

2i  va/m>=Ym  * 

Y lLN*2)*Yf I J+XC  * 


IF  { JIK.6Q.1.AN0.J8.EQ.1)  GO  TO  40 
NV1*NY~£ 

Y I LN*3 ) -Y I  NY ) ♦( Y (NVi )-Y( NV ) ) »XC/W 
GO  TO  4i 

40  W i^TWGK 

mN+3I  =  {W+Wl)*{Y(NV)-Y(NV-lH/Wl+YCNV-i) 

41  YfLN+4)=i,-{Y{LN+3)/V0)**2 
LN=*LN+§ 

90  V=0* 

73  CONTINUE  * 

58  IF  { LN*GT « 17983)  GO  TO  100 

RETURN  -  .  2 

100  WRITE  (6,601) 

STOP  . 

END  .... 

83 
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$09 ROUT INS  SU853  (Vf>tNt 
DIMENSION  f U  l),YU8000)fSLxNH3C30) 

COMMON  /COM/  UTtWl 
COMMON  /PLOT/  RANK 

COMMON  BtNXC,NYR:>NKXVC,NKYV€,NDXC*NDYCjNCFCT,NBETACtKJ*8IT$f ISTt  1 

IHONEiHTMO.ONEK, TWOK,If IN* J»JN , !XJ,XI J, JiK, I JK, IT, NDX, NDY,NFUNCT,Y9  l 

2SlIN£S,HI,H2,TKl,TK2,NBI»N82fNB3,N64»NCI#NC2fNC3,NC4»Nl,N2tU0,Ul,  1 

3U3  »U2«U4» |PRINiNRS0U£#NSL»NSLNC»NGRDC»N5RDXC*NGRDYCt A  » YH»NGRAD»  I 

4NGRADX* NGRADYtLNj CSCH»NG*  EXTRAP, SSRtK,  NSWPS, OMEGA, CVGS, CONS? »  l 

5STReAMtKNrR»A,C,D,E,FfG,HKl,HKl2iHI2sHK2*8Y,BX9NKXV,NKYV,N3.  i 

6P,0tR,S,T,V,W,XC,YC* IH  1 

YIINJ-SUmgSINI  1 

Yf INI -RANK 

P«$UN£SIN)-U6  l 

S-.UO-Sl  *NE$(N J  1 


IP  {A8$!P*Q'SlT*ASS{G-?n  GQ  TO  2 
P-SUNESfNJ 

IF  ' AXYH.€Q.O.O.OR*$%£0« U I  60  TO  20 

XC- ( P-UO ) * { P-U6 )/{{ U7-U0 ) » i U7-U6  H  *  Wl ♦ { P-U7 ) ♦ ! P-UO) / { ( U6-U7 ) *  f  U6 


IF  { A8S i XC ) « LT«  A3S I W  J )  GO  TO  21 

20  XC*  I P ~U0 )  /  * U6-U0.I 

21  V*i*  . 

2  RETURN. 

END 


1 

1 

1922 
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SI8FTC  WAREGN 

SUBROUTINE  REGION  0. 

DIMENSION  8m,Y!  18000)  tSLlNESI 30) 

COMMON  B,NXC,NYR,NKXVC,NKYVC,NDXC,NDYC»NCFCT,N8£TAC,KJ, BITS, 1ST,  O’ 

1H0NE , HTWO ,ON£K, TWOK , I , IN t J , JN , IK J ,K IJ , J I K, I JK, I T , NDX, NDY, NFUNCT , Y ,  O' 

2SLINES,H1,H2,TKI,TK2,NB1 }NB2,NB3,NB4,NC1,NC2,NC3,NC4,N1,NL,U0,U1»  0. 

3U3,U2,U4, IPRIN,NRSDUE,N$L,NSLNC»NGRDC,NGRDXC,NGRDYC,AXYM,NCRAD,  0. 

4NGRA0X,NGRADY,LN,CSCH,NG, I XTRAP,SSRCH, NSWPS, OMEGA, CVGS, CONST y  0 

5STREAM, KNTR, A*C»  D»  E»  F»G»  HK 1 » HK12»H12,HK2» BY, BX  »NKXV»  NKYV»N3,  0.' 

6P,Q,R,S,T,V,W,XC,YC, IH,NTAPE  0: 

COMMON  U8DY,LB0Y  0: 

226  FORMAT  { 3H  A=F12.6/3H  C=F12.6/3H  E=F12.6/4H  Y1=F!1*6/4H  Y2*F11.6/5 

X3H  THE  DIFFERENTIAL  EQUATION  IS  NOT  ELLIPTIC  AT  IY1»Y2)/57H  IF  SUB 
^ROUTINE  ACDEFG  IS  USED.  OTHERWISE  CHECK  A,C,AND  F.//) 

C  STATEMENT  112  TO  4  SET  MESH  VALUES  FOR  DERIBATIVE  ENFORCEMENT  0: 

IF  (  IXTRAP )  112,113,112  0.' 

112  I KJ=0  o: 

kij=o  o: 

j*ik=0'  o: 

i  jk=o  o: 

HON£=Y ( NDX )  0; 

HTWO=Y  ( NDX- 1 )  0! 

ONEK=Y I NDY )  01 

T WOK~Y ( NDY- 1 )  0! 

IF  ( UBDY }  1,2,1  0: 

1  ONEK=TWOK  Oj 

GO  TO  3  0| 

2  IF  ( L8DY )  4,113,4  0! 

4  TWOK=ONEK  0! 

go  to  3  o: 

113  CALL  GRDSTL  0! 

3  IF  (CONST)  141,149,141  0) 

C  CONST,  AN  INPUT  VALUE,  CONTROLS  USAGE  OF  ACDEFG.  0! 

C  XC,YC  ARE  X  AND  Y  COORDINATE  VALUES.  0; 

149  CALL  ACDEFG  Oi 

141  IF  ( A  *  C )  147,147,43  0! 

C  STATEMENT  141  CHECKS  ELLIPTICITY  0! 

147  WRITE  (3,226)  A, C, £, Y( I ) , Y< J I 

KNTR  =  l  0: 

43  CALL  CNSTNT  0; 

C  SUBROUTINE  CNSTNT  CALCULATES  COEFFICIENTS  OF  DIFFERENCE  ANALOGUE  0) 

75  RETURN  0; 

END 
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SIBFTC  WnSORT 

SUBROUTINE  SORT (KASE) 

0 1 PENS  ION  XTI250) ,YTI 250) 

DIMENSION  XP(250),YPI250J 
01  MENS  ION  TAI4G00) 

DIMENSION  TPI91 

DIMENSION  B{ 1) ,Y  0.8000) ,SLINESI30) 
DIMENSION  ITLEI24),HEAD(2».3) 


COMMON  B»NXC».NYR».NKa'.'C*.NKYVC f.NDXC»NDYC»NCFCT »NBETAC»KJ»BIT$»IST*  1 

1H0NE » HTWO*  ONEK »  TWOK , UI N » J , JN , I K J, K I J , J I K ,.1  JK,.I T » NDX , NDY.NFUNCT , Y »  i 

2SLINES,HI,H2,TK  WTK2„NBl,Na2,NB3,.NB4,NCl,NC2,NC3,NC4,Nl,N2,U0,Ul»  l 

3U3,U2,U4rIPR  tN,-NRSOUE,.NSL,.NSLNC  »NGRDC  »NGRDXC  ».N6RDYC , AXYM, NGRAD,  i 

4NGRACX,.NGRADY,LN,.CSCH,Nf?,NDER  I  VtSSRCH,  NS  WPS,  OMEGA,  CVGS,  CONST  *  1 

5STREAMf  KNTR ,  At C »  Dt-E »  Ff-G*  HKil , HKL2*.H12»HK2»3Y t BXfNKXV»NKYV|N3»  l 

6P,Q,R»9tT,V,W,XC,YC,  IH#.NTAPE  l 


COMMON/ PRSCQF/  FIRST, VO 

COMMON  /C0M2/  XSCL ,  YSCL ,  PUNCH,  I  TIE  ,  YO,  YUXl  1 ,  XL2  ►YU  ,  YL2  *XSCLA, 
IYSCLA,XLIA»YLlA,XL2AfYL2A 
COMMON  /C0M3/  YMAX 
EQUIVALENCE  ( T  A ( I ) ,Y I 1 1000) ) 

EQUIVALENCE  (  XT  (  I )  ,.XP(  1 ) ) , ( YT( 1 )  ,YP  III) 

600  FORMAT  ( 1H W2A6 , 1H  =  F 14 . 7//8X , 1HX,-1 3X , IHY ».10X»>8HVEL0C I  TY,6X# 
l8HVEL0CITY,6Xf 8HPRESSUR E/47X, 5HR AT  10, 9X , 11HC06FF IC I  ENT/ ) 

620  FORMAT  (IX,5F14.2) 

630  FORMAT  (1X^131 
640  FORMAT  IIX//L 
700  FORMAT  (5E14.81 
1001  FORMAT  ( 5F 14 . 7 1 

DATA  C  HE  AO  (  1 ,11,  1=1,2)/ 6H  PSI,6H  / ,.( HEAD  (  U2  )  ,  I  *1 ,2  ) /6HX-06 

XRI ,6HVAT I VE/,1  HEAD! U3) , 1  =  1 , 2J /6HY-DERI , 6HVATI VE/ 

R ANK= i • 

LN=LN-5 
XL IP=XL1 
XL2P=XL2 
YL 1 P=YL 1 
YL2P=YL2 

IF  (XSCLA.EQ.O.)  IST0P*1 

L INE=0 

TP11H23. 

TP { 2 ) fc4« 

TP  C  3 ) =XSCL 
TPI4HYSCL 
TP ( 5 ) k0« 

TP ( 6 ) fc0« 

TP  I  8 ) fc4» 
rP(9)«=0. 

2  NSTOR  =LN*5 
•I  STOPfeO 

5  00  20  I=NGRDG,.LN,5 

IF  (Y! I ) .NE.RANXi  GO  TO  20 
.10  Y I NSTOR  )  SYI  U 

Y(N$TOR+l)nY(  HI) 
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Y (NSTOR+2 ) { ( ♦2 ) 

Y(NST0R+3J-YUt:3l 
Y  ( NSTQR+4) ^YC 1+4) 

NST0R*NST0R+5 
20  CONTINUE 
GO  TO  23 

22  RANK=RANK-1,0 

23  NSTORfcNSTQR-5 
NSTCR.l*LN+5 

*F  ( NSTGR-NS  TOR  I  *LE . 6 )  GO  TO  102 
46  LCNE-O* 

XL 1P=XLI 
XL2P=XL2 
YLIP’YH 
YL2P=YL2 

971  WR3TE  (6,600.)  (  HcADU  tKASE  ) ,  1*1, 2) ,  SLINES(NG) 

DO  775  I=tN$T0RL,NST0R»5 

IF  (Y(t+1).EQ.Y(  1+6)  *AND.Y{  l*.2)»EQ»  YII+7))  GO  TO  775 
>VR=Y(Hi3)/V0 

WRITE  (6,1001)  VM  +  l),Y(I+2),YU+:3),VR,YU.+4) 

LINE=LIN£*1 

<IF  (  LINE  *LT. 50L  GO  TO  775 
LINE=Q 

WRLTE  (6,600)  ( HEAD i J ,KAS6 ) , J ~1 , 2 i , SU NES ( NG ) 

775  CONTINUE 
72  J*0 

DO  100  I -tNS  TORI,  NS  TOR  ,5 
JfcJ  +  1 

XPIJ)tY(IH) 

100  YP( J)*Y( !+2) 

KNTR  =0 

00  65  JK*.ltJ 

IF  (XP( JK) .LT.XLlP.OK*X<M JK ) * GT. XL2P. OR. YP ( JK) .IT. YLLP.QR. YP( JK)  . 
XGT.YL2P)  GO  TO  62 
GO  TO  65 
62  XP(JK)~BIT9 
KNTR=KNTR*1 

65  CONI' NUE 
KT  =  0 

CO  68  JK=?  1 » J 

IF  (XP( JK) .NE.BITS)  GO  TG  68 

66  KT=KT* 1 

jkt=jk+:kt 

IF  ( XP ( JKT ) . EQ. B  ITS )  GO  TO  66 
DO  67  L  J^JKi,  J 
JKT=l J+KT 
XP(LJ)=XPl JKTI 

67  YP(LJ)=YP(  JKH 
KT  30 

68  CONTINUE 
JfcJ-KNTR 
AJoJ 
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TP  I  7} *=AJ 

IF  IISTOP.NE.Oi  GO  TQ  80 
WRITE  (8)  .1 

WRITE  (8)  I XP 1 1 )  ►  I  =  1  »J ). »  ( YP { I )  »  i;^t »  J ) 

XL  iP  =  XL IA 

XL2P=XL2A 

YL1P=YLIA 

YL2P=YL2A 

I STGP*2 

GO  TO  72 

80  IF  USTOP.EQai  GO  TO  81 
TP  ( 3 ) =XSCLA 

TP  ( A  )  =YSCLA 
TP ( 2 ) 

81  CALL  GPL0T9!TP,.XP,.YP) 

RANK®  RANK*1. 

GO  TO  2 

102  i F  (NG.LT.NSU  RETURN 
TP (  7 )  *=0 
TP ( 9 ) =-l • 

TP( 2 ) -4* 

CALL  GPLOTS  !TP,XP,YP) 

IF  (-I5T0P.NE.2L  GO  TO  120 
REWINO  8 
TP  { 3  )  =XSCL 
TP  ( 4  )  =YSCL 
TP( 9 ) fc0. 

DO  101  I ~ 1 » NSL 
READ  (8)  j 

READ  (8)  (XP(ML,M  =  l,J),(YPlM).'sM*i:f  J) 
TP  (  7  )  =  J 
TPt  2) =4. 

101  CALL  GPLOTS  (TP,XP»YP)- 
TP(  7)5=0. 

TP(  9)  at-1. 

TP  t 2)®4. 

CALL  GPL0T3  ( TP »  XP » YP ) 

120  CALL  REMOVE! 23*3) 

IF  (PUNCH. EQ.O. )  GO  TO  99 
YT (  l )  =Y( IST+i) 

CO  123  I  *2 » NYR 
ID=NKYV+I-1 

123  YT!  I  )=Y!  ID)+YTU-1): 

XT!  1 )  *Y (K J  ) 

00  124  I  =2 » NXC 
ID=NKXVf  I~1 

124  XT!  I  }=Y(  lOi+XTU-U 
t a ( n*i. 

TA( 2 ) kNYRf l 
TA(3)*=NXC  +  l 
T  A  ( 4  )  **0. 

DO  125  I  3  l » NYR: 


D3~696l^ 


125  TA( I f 4 ) =YT { J ) 

NFUNCr=NCFCT-.l 
JK=NYR+:5 
DO  256  1*1, NXC 
TAf JK)=xr(l) 

JK= JK+ 1 

00  255  NQ*l,NYft 
NFUNCT=.NFUNCTH 
TA l JK ) =Y( NFUNCT) 

JK«JK+1 

255  CONTINUE 

256  CONTINUE 
JK-JK-l 

WRITE  17,700) (TA(I) , Ini, 3) 
WRITE  (7,700)  (TA( I ) ,I=4,.JK) 
99  RETURN 
END 


) 
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I 


HBFTC  WAMARK 

•MARC  1 

SUBROUTINE  MARC  I 

*  THIS  SUBROUTINE  PERFORMS  THE  FUNCTION  OF  SEARCH  IN  LINK  3.  I 

D I  PENS  I  ON  B  ( 1 ) » Y { 18000) , SLINESI 30) 

COMMON  8,NXC,NYR,NKXVC,NKYVC,NDXC»NDYC,NCFCT»NBETAC,KJ,BITS,IST,  l 

lHONE,HTKO,ONEK,ThOK,  I , IN, J, JN, IK J , K I J , J IK , I JK , I T , NDX , NDY , NFUNCT, Y,  l 
2SLINES,H1,H2,TKI,TK2,NBI,N82,NB3,NB4,NC1,NC2,NC3,NC4,NI,N2,U0,U1,  I 

3U3*U2»U4»IPRIN,NRSDU  £ ,NSL » NSLNC , NGRDC , NGRDXC  »NGRDYC , AXYM,NGRAD»  l 

4NGRADX , NGRADY, LN, CSCH,NG, EXTRAP ,  S  SR  CH  (>  NS  WPS ,  OMEGA,  CVGS,  CONST  ,  1 

5S fREAM.  ,KNTR,A,C,0,E, F»G,HK1,HKI2,HI2,HK2»BY*BX»NKXV»NKYV*N3*  1 

6P,Q,R,S»T,V»W,XC»YC» IH*JH,KH  l 


COMMON  /VEL/  NGROVC, AIR 
COMMON/BDRY/  JB 
COMMON  /ZCOM/  Z 
COMMON  /LNK5/  NV 
COMMON  /PLOT/  RANK 
600  FORMAT  { 1 12 ) 

620  FORMAT  (1X,6F14.2) 

*  INITIALIZE  ALL  SUBSCRIPTS  TO  PICK  UP  INFORMATION  IN  Y  ARRAY.  I 

•  1 
HL=Q.O 

SRANK* 1.0 

FLAG=0.0 

Z=0. 

LN=NGRDC 


NGRACX^NGROXC-l  l1 

NGRADY=NGROYC-l  l 

NV=NGROVC-I 

NKXV=NKXVC-l  I 

NFUNCT =NCFCT-1  l 

NDX=NDXC-l  1 

J  =  KJ  l 

*  l 

*  SET  LOOP  TO  PROGRESS  FROM  COLUMN  TO  COLUMN.  I 

*  l 
FST=0«. 

DO  76  NX=I ,NXC  l 

*  CHANGE  SUBSCRIPTS  AFFECTED  BY  A  CHANGE  OF  COLUMN.  I 

RANK=SRANK 

*  l 

NDY=NDYC~T  l 

NKYV’NKYVC-l  l 

NDX=NDX*l  l 

I  =  ISTfl  l 

NKXV  =  NKXVU  I 

*  1 

*  SET  LOOP  TO  PROGRESS  ALONG  A  COLUMN.  I 

DO  75  NY  =  I » NYR 

*  CHANGE  SUBSCRIPTS  AFFECTEO  BY  A  CHANGE  ALONG  A  COLUMN  1 

*  l 

I LN-LN 
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NGRACX-NGRAOX+1 
NGRACY=NCRADY+l 
NFUNCT  =NFUNCT  +  l 
NV=NV*1 
NKYV=NKYV+i 
NDY=NDY+ I 

SET  SUBSCRIPTS  TO  FIND  BOUNDARY  INFORMATION  FROM  THE  INPUT  ARRAY, 

JN=*  J*2 
IN=I+2 

KYVA=Y(NKYV)/4. 

KYVB  KYVA 

LCCP  10  AND  112  ARE  NECESSARY  FOR  COMPLEX  BOUNDARIES. 

DO  10  IS=l,KYVA 

LOCATE  POSITION  OF  POINT  IN  THE  ENTIRE  SYSTEM. 

THE  POINT  WILL  EITHER  86  INTERIOR, EXTSR IOR  ORBOUNDARY. 

IF ( Y ( IN)~Y{ J) ) 1,80,75 
1  IFCYI  IN  +  2)-YU) )  110,80,7 
7  KXVA=Y(NKXV)/4. 

KXV8=KXVA 
DO  112  JS= 1 »  KXVA 
IF(YCJN) —Y { I ) ) 11,75,75 

11  IF(Y(JN+2)-Y( I) ) 12, 75,13 

12  KXV8-KXVB- l 

IF  (KXVB)  75,75,112 
112  JN=JN*4 
llO  KYVB  =  KYVB- 1 

IF  (KYVB)  75,75,10 
10  IN- IN+4 

13  CALL  GRDST2 

14  CALL  SGSRCH 

IF  (LN.GT.ILN)  RANK=RANK+1 « 
l- 0. 

75  I  =  I  +  I  NT ( Y  C NKYV ) ) +  2 

IF  (RANK.LE.2.0.0R.FST.NE.O.)  GO  TO  50 
IF  (FLAG. EQ. 1,0)  GO  TO  60 
IF  (AIR. NE. 1.0)  GQ  TO  55 
52  FST=1,0 
HL  =  l  •  0 

IKK=5.0* ( 2.0+ FLAG) 

I OX=FLOAT ( LN ) - ( 10.0+ FLAG* 10.0) 

00  57  I K  = 1 , IKK 
K0NT=LN+6-IK 
57  Y(K0NT)=Y(K0NT-5) 

Y  I  I  OX  )  =*2.0 
GO  TO  50 
55  FLAG3 1.0 


GO  ro  25 
6C  CONTINUE 

IF  {RANK«EQo3* )  GO  TO  52 

SRANK==2.0 

RANK=SRANK 

FST=l.O 

IF  ( AIR.LT .1,0)  GO  TO  74 
DO  65  IK" i» 20 
K0NT=LN+6-lK 
65  Y  (  KONT  )  =Y.{  KONT-5  ) 

Y ( LN-30 )  =2.0 
Y ( LN-25 ) =3.0 
Y(LN-20)=i.0 
Y ( LN- 15 ) =4.0 
Y { LN- 10 ) =2.0 
Y (LN-5 ) =3.0 
Y I LN ) =4. 0 
LN=LN+5 
GO  TO  25 

74  DO  745  I K= 1 1  25 
K0NT=LN+6-IK 
745  Y { KONT ) =Y { KONT-5  ) 

Y ( LN-25  )  =2.0 
Y ( LN-20 ) *3*0 
Y (LN~15)=4oO 
Y(LN-t0)=2.0 
Y ( LN-5 ) =3.0 
Y { LN ) =4. 0 
LN=LN+5 
GO  TO  25 
50  CONTINUE 

If  (RANK. NE. 4.0)  GO  TO  25 

IF  (HL.6Q.1.0)  GO  TO  25 

SRANK  =  5 . 0 

RANK=SRANK 

DO  24  IK  =  1,  10 

KONT-LN+ll-IK 

24  Y ( KONT ) =Y ( KONT-IO) 

LN=LN> 10 

Y ( LN-10 )  =R-*NK 
Y(LN-5)-RANK+1.0 

25  CONTINUE 

76  J  =  J  + 1 N  f  ( Y  {  NK  X  V )  )  ■*■  2 

26  RETURN 
80  2=1. 

GO  TO  7 
ENO 
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! 

i  • 


U&ti'i/  X)Asi}£e,t~  Cf'ULj&etotei 


REVLTR: 


SECT 


8-30*3  R f 


NO. 

REVLTR: 

SECT 

PAGE 
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e*jos>  r 8 


REV  LTR: 


my 


(jgAb  mus} 


(SR 


/Q\ 

VI 71 ' 


(ft£AD  //BAD 


BLANK 


) 


hot  BLANK 


read  PICT 
°ARAmr£RS 
] 


D 


Read  alwc.p .w£\ 

UPlNF.  PA.  WA  ) 

. . 


HMD  XIWtfftM  .Xi)J?sA 
J>£LYDP  XL  m  i  J 


vy/kz/Vj  />  pgcr,  flwmt  j 


^EAD  m/fi/LjAUMDA^S 
\RE3IL  mu.  LAm\R  ) 

('mJdnAmi  D&wt\ 

\&£M ,  /.?£,*/.  XPSFZ  J 

- /- . . — ^ 


sect 
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PAGE  1*56 


E-J033  R  I 


TL3  -  .Vrj 
7Xy  =  .7/ry 
ris  *  -s/ts 
TLi  -  Ti-.sAi 
Til  -  TsfZ.Ufc 


'ALPHA  ( TLHEAD 


Tt-0 


tmois  v 

(t,xplt.mj) 


I  NO.  P3-09Q. 
[sect  PAGE  157 


REVLTR: 


■  *-  »-•  •  *«  ■  I  • 


:  *  •  i 


) 


5 


) 


{RE AX>  mAPI  YTMPz\ 
'\XTRAPi,  XTfiAPZ,  TABJLEJ 


TA81 6)  -=io 
TASKS)  *Z.O 
TASKS)  -JOO'O 
TASK'D  *  & 
TA3lCs)~0. 
HK^S 


NK-HSt! 
m^mti 
TAB  I hh  -  TAB6>fj^  ’''AML 
NK-  NKt  I 
a  NKl  A  l 


a 


=  0 


TAB  I ~  o.  0 


'3 


NO. 
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SECT 

PAGE 

158 

e-Jon  ni 


/&£7&*A/£?\W0.  W-696I 

REVLTR:  ■ 

SECT  JPAGE  160 

E-JOJJ  Rt 


NO. 
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SECT 

i 

PAGE 
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REVLTR: 


▼  * 


1 


•  •  i 


r  l 


r 


e-joij  r t 


&£7JZrjr/y£Z 

NO.  D3-®61 

REVLTR: 

SECT 

PAGE  166 

REVLTR: 


J-J-hNYR  ;i 


[71 C  TA33jXiCCOO)  +  0.S 
I*1  10000 


ThBh  -  X2 


I-I+l 


>NXC, 


kYMAX 
TeLK  =D£lXn 


m\£^£.££_  oax-  '>Ein 


tRbJFl 


>  XREFZ 


REVLTR: 


e-sojs  ns 


D3-^96l 

[sect  '  jpAGTTi™ 


SECT 


*  i 


■v 


£3 

NO.  D3-6961 

4: 

') 

REVLTR: 
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PAGE  173 

COSH  «S 


f TLUMS<jrA&i\ 
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'tluu$  C 7am\ 
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1 
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>0 

(248) 

1 

, _ 1 _ 

<T/VP/  y 

'y 

5** 

im 


PHIU  -  0 


'tIjUYYS  (EAQSj  \ 
UljjULjC/HItt) 


REVLTR: 

e-ioji  «i 


>0  X  >v 

— — <  IUD%  > - 


phluL  =  / 


Z>V£2, 


/UUXYSirABJj 

m,y uhtotPHxoC) 


I pwr-o  ' 


NO,  I 

SECT  PAGE 
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VBV1L  =(UL£ui)t>em 
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'deivr 


\PmR=i.o\ 


JND3 


IND3 


\PHllR-0 
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ind y 


JyYD  V 
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'nwi  (rA33.  N 
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— 

m _ - 
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vw?=  WWR-mLA 
p^y/? 


3(o 

VH  =  0  j 

J 

28S 


Tiwcx^  (M55, 
tlLYLLO^mi) 

v  '  .  '  v 


)0 


NO. 
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SECT 

PAGE 
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e-iosj  «i 


REVLTR: 


£7£?jZrS/V£?  1  NO.  D3-6961 
SECT  PAGE 181 


C-JOJJ  At 


E«»033  n  1 
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T 

>0  ^  y 

ZONTER^dORTER-ii 

X 
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>0 
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M 

A 
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NO.  D3-6961 

SECT 
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T 


±.0! 
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V 
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- < 
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3011  \ 

/ 

ho 
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HO 
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DvytA  =  VDYWMI  -  VPyivW 


NO.  D3-6961 

REVLTR: 

[SECT  j  PAGE  187  ~ 

E-301J  R  t 


>  <5 


vayww=vryvfAi 


Dvy/w-i>vy/\ 
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i>^r _ 
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m 
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n 
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10.4  Appendix  4  -  Water  Droplet  Trajectory  Computer  Program  Listing 
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NC.  D3-6961 

REVLTR: 

SECT 

PAGE  206 

E-3033  R  I 


o  o  o 


.  >• 


SIBFTC  DRCP 

DIMENSION  T{7) #XPL0T(5  50) ♦YPLOTC  550) tXPLT(3lOO) ,YPLT( 3100) 
DIMENSION  TAB1( 100) »TAB2( 100) »TA33(  5C00) , TAB4 { 200> * TA85 ( 200) 
DIMENSION  TAB6(200) ,7AB7{200) 

DIMENSION  HEAD! 12) # TA{ 7) , IPCt 100) 

DIMENSION  TL(9) 

EXTERNAL  GOTHIC 
COMMON  TAB3 

DIMENSION  TAB6I200) ,TAB7(200) 

DIMENSION  T(7) ,XPL0T(200) , YPL0T(?00) ,XPLTI20G0) ,YPLTt2000> 
DIMENSION  HEAD! 12),TA(7) , I  PC  C 100) 

DIMENSION  TAB3I 100) ,TAB9{ 100) 

INTEGER  HEAD 

DIMENSION  TAB 10 ( 500 )» TAB 111  800) 

DIMENSION  TL(9) 

REAL  M.Q,MOT,  LANGMR 
C  COMMON  TAB3 

C  GOTHIC 

10CC  FORMAT! 6E 12. 6) 

8787  F0RMAT(4X,4HDTW=El2.4t4X,5HV0YW=E12.4»4X*7HV0YWNl=E12.4/4X,5HVDXW= 
1E12.4,4X, 7HV0XWN1=E 12.4, 4X,3HAK  =  E1 2.4/4X,3HCK  =  E«.  2.4, 4X,4HREN  =  E1 2.4 
2,4X,6HREINF=E12.4/4X,4HVDX=E12.4,4X,4HVDY=E12.4,4X,3HVX=E12.4/4X, 
33HVY=E12.4) 

1240  F0RMAT(4Xf7HWATRAT=E20.8,7X,4HAVS=E20.8/4X,3HXP=E20.8,llX,3HYP=E20 
1 . 8  »  9X  ,  5HWATER  ,F14.6) 

DIMENSION  Y(  150) ,X(  150) 

6C0  FORMAT  (6E1S.5) 

601  FORMAT  (IX) 

603  FORMAT  {  lXt 4F14.5) 

604  FORMAT  t /4H  XP  =  F1?.6/5H  XN1=F11.6/) 

610  FORMAT  ( 9H  AHLSTEDT ) 

63C  FORMAT  (  1 X  ,  2HR  =  F  1  2 * 6 / / 6 X  ,  3HAVS  ,  1 IX  ,  6HWATRAT/ /  (  F  1 4 . 6  » 2 X  ,  F  6  )  ) 

635  FORMAT  (///6X*3HAVS, 1 1 X , 6HCUMWA T // ( F 1 4 . 6 , 2X , F 1 4. 6 ) ) 

640  FORMAT  {  1X,6F14.6 ) 

65C  FORMAT  (  1X,8F9. 3) 

660  FORMAT  (17H  PARTICLE  TRAPPED) 

661  FORMAT  ( 38H  TOO  MANY  POINTS  IN  PLOT  OF  WHOLE  AREA) 

662  FORMAT  (40H  TOO  MANY  POINTS  IN  RESTRICTED  AREA  PLOT) 

670  FORMAT  ( 8H  l  T  ABLE  I/33H  X=OISTANC£  TO  RIGHT  OF  H I GHL I GH T/45H  Y  =  LOWE 
XR  CCWL  SURFACE  OISTANCE  FROM  H IGHL I GHT// 1  OX , 1HX , 2  OX , l HY/ ) 

671  FORMAT  '£H  TABLE  2/33H  X=DISTANCE  TO  RIGHT  OF  HI GHL I GHT/46H  Y=UPPE 
XR  COWL  SURFACE  DISTANCE  FROM  H IGHL I GHT / / 10X , IHX , 20X , IHY/ ) 

6C5  FORMAT (IX//// ) 

673  FORMAT  ( 8 H  TABLE  4/34H  X=Y  VALUE  OF  POTENTIAL  FLOW  F I ELD/53H  v  =  FIR 
XST  BOUNDARY  VALUE  OF  X  AFTER  LEFT  HAND  BOUNDAR Y/ / l OX , 1HX ♦ 20X , 1 HY/ ) 

674  FORMAT  (8H  TABLE  5/34H  X=Y  VALUE  Or  POTENTIAL  FLOW  FI6LD/54H  Y=SEC 
X OND  BCUNCARY  VALUE  OF  X  AFTER  LEFT  HAND  BOUNDARY/ / 1  OX ,  1HX , 20X , IHY ) 

675  FORMAT  (8H  TABLE  6/34H  X  =  X  VALUE  OF  POTENTIAL  FLOW  FISLD/50H  Y  =  F I R 
XST  BOUNDARY  VALUE  OF  Y  ABOVE  LOWER  BOUNDARY/ / l OX ! HX , 20X , IHY/ ) 

676  FORMAT  ( 8H  TABLE  7/34H  X  =  X  VALUE  OF  POTENTIAL  FLOW  r I E  LD/50H  Y  =  SEC 
XOND  BOUNDARY  VALUE  OF  Y  ABOVE  LOWER  BOUNDARY// 1  OX  1 HX , 20X 1 HY/ ) 

677  FORMAT  (8H  TABLE  8/34H  X  =  X  VALUE  OF  POTENTIAL  FLOW  F I  ELD/4 1H  CENTE 
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XRLINE  OR  LOWEST  BOUNDARY  VALUE  OF  Y. / 10X , 1HX , 20X , IHY) 
o78  FORMAT  { 3H  T ABLE  9/39H  X=DISTANCE  TO  RIGHT  OF  HIGHLIGHT  (XHLJ/63H 

1  Y  =  CEnT£R8GDY  SURFACE  DISTANCE  FROM  CENTER80DY  HiGHLIGHT  { XHL l ) // LO 

2  X , 1HX ,  20X  » IHY ) 

6fiC  FORMAT  ( 10X,F10.5, L0X,F10.5) 

3  ICO  FORvAT { 1 2 A6 ) 

3101  FORMAT { 1HI, 12A6) 

3102  FORMAT ( IhO, 10X4hALWC17XlHR14X4HtJINFI3X5HUPINF16X2HPA16X2HWAi 

3104  FORMAT  (IHO, l 1 X , 3HYSL , 1 5X, 3HYSM, I 3X , 5HVXWI N , 13X , 5HVYWI N , 13X , 
16FDIRECT  r 13X »  6HRAND0M  ) 

3105  FORMAT  { 1HC,  l IX , 3FXHL  » 1 5X,  3HYHL  , 12X , 6HALAMDA , 1 3X, 5HRES I L , I4X , 4HRH0. 
XW,  12X,6HLANGMR) 

3165  FORMAT  ( 1H0,9X, 5HREINF, 16X,2HAK) 

3IC6  FORMAT  ( IX , 6 ( E 1 5 . 7 , 3X )) 

3107  FORMAT! IHO) 

31C8  FORMAT  ( IHO, 10X,4HXMAX, I4X , 4HYMA X , 1 3 X , 5H YDROP , 1 2X , 6HDELYDP , I6X  , 
12HXL , 1 4X , 4HXHL I  ) 

31C5  FORMAT  (2X,4HXNNl ,5X,4HYMN1 , 6X , 3H0T W , 6X , 2HCK , 7X , 3HREN , 9X , 3HVDX , 

17X, 6HVDXWNI , 6X , 6HVDYWN I , 8X , 3HVDY , 6X , 3HXKT , 3X , 6HD I REC T ) 

3 1  I C  FORMAT  (IH  »F7.3,Fl0.4,lX,F8.5»iX,E8.1»lX,E11.3»lX,E11.3»lX,£ll.3» 
IIX,E11.3»1X,E11.3»1X»3F5.1) 

3111  FCRmaT(8H0V0YWM  =E l 2 . 4 , 5X 7HV0XWM  =£12.4) 

3112  FORMAT  (6E12.6) 

3113  FORMAT  (5E14.6) 

3114  FORvAT  { 1H0,8X,6HGRAVTY, 13X, 5HDELXI , 1 3 X , 5HDE LX2 , 1 3X , 5HXREF 1 , 1 3X , 
15HXREF2) 

3115  FORMAT  ( 1H0,8X,6HYTRAP1, 12X , 6HYTRAP2 , 1 2X , 6HXTRAP l , 1 2X , 6HXTRAP2 , l 2X 
X,5HTAGLE) 

DATA  IBLANK/6H  / 

L ANGMR=0. 

REw INC  23 
WRITE  (23,610) 

10  =  0 

CO  2009  1=1,800,2 
FILL  =  FLOAT(  1-1)/*. -30. 

2  C  C  ?  TAQ  11(1) =F  ILL 
BCU\CE=0 .0 

2001  CALL  RENTRY  (N) 

IF  (N.NE.O)  GO  TO  6041 
2CCC  CONTINUE 

C  ^L  L  REMOVE  (23,3) 

TL ( 2 )  =  0.0 
T  L ( S )  =  10.0 
T  L ( 9 )  =  -1.0 
GO  TO  2004 
6041  X55=X55 
2003  T  A ( 7  )  =  10 

I  F  (  I C  -  2  )  3394,3394,3393 

3393  CALL  MPtOTS ( T A , XPLO T , YPLOT  ) 

3394  T A ( 7 )  =  0.0 

T  L I  1)  =  T A  (  1) 

T  L ( 3 ) =  . 3/T  A ( 3 ) 
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TL(4)  *  .7/ TA  {  4  ) 

1"L(5)  =  *  5/ T  A ( 3  ) 

U(6}=  TAC6)-,5/TAC-45 
TL( 7}  =  7 A C  5 }  +  2 . 25/ TA { 3 ) 

CALL  ALPHA! TL, HEAD, 72, GOTHIC) 

CALL  .vPLOT$  {  TA,  XPLOT,  YPLOT) 

1(1)  =  TA { 1 ) 

T  12)  =  TA ( 2  J 
IF(JP)3396, 3396, 3395 

3395  IP  s  ip  +  j 
1PCIIP)  =  JP 

3396  IF( IP)2002, 2002, 3397 

3397  IPA  =  1 

f)0  3398  1  =  1, IP 
J  =  IPCU) 

T  {  7  )  *  J 

CALL  KPLOTS (T,XPLT( IPA),YPLT( IPA) ) 

3398  IPA  =  IPA  +  J 
TL(3)  =  . 3/ T ( 3 ) 

TL(4)  =  .7 /I(4) 

TL(5)  =  .5/TI3) 

TL(6)  =  T ( 6 )  -  . 5/T { A ) 

TL(7J  =  T ( 5 )  +  2.25/TC3) 

CALL  ALPHA{ TL, HEAD, 72, GOTHIC) 

T  (  7  )  =  0 

CALL  .MPLOTSI  T,XPLT,  YPLT) 

IF  (LANGMRoGT.O. )  GO  TO  2005 
GO  TO  2002 

20CA  READ  (5,3112)  ( TA81 ( J ) , J=1 , 3 ) 

I  A  =  T  AB 1 ( l)*TAQl(2)*rABlf3)+3.0 
REAC  (5,3112)  (TAR1( J),J=4,  IA) 

REAC  (5,3112)(TAB2(J),J=l,3) 

I A=  TAB2 ( 1 ) »r  AB2( 2) *FAB2( 3) +3.0 
REAC  (5,3112)  ( TA82 ( J ) , J=4, IA ) 

REAC  (5, 3113) (TAB3( J } ,J=1,3) 

I  A=  r  AB3 {  1  )*TAB3( 2 ) * TAB3 ( 3) +3. 0 
REAC  (5,3113)  (TAB3(J),J=4,IA) 

RF AC  (5,3112) (TAB4(J),J  =  lf 3) 

I A  =  FAB4 ( i)»TAB4(2)*TAB4(3)+3.0 
REAC  (5,3112)  !  TAB4{  J),.J  =  4,  IA) 

REAC  (5,3112)  (TAB5(J),J  =  1, 3) 

I A  =  FAB5 (  l)*TA85(2)*TAB5I3)+3.0 
REAC  (5,3112)  ( T AG5 ;  '), J=4, IA) 

REAC  (5,3112)(TAB6(J),J=!,3) 

IA  —  TAB6( 1)*TAB6(2)*TAB6( 3) *3.0 
REAC  (5,3112)  ( TAB6 ( J), J=4, IA) 

REAC  (5,3112)(TAB7(J),J=1,3) 
IA=FA87(l)*FAB7(2)»TAB7(3)f3.0 
READ  (5,3112){TAB7{J),J=4, IA) 

REAC  (5,3112) (TAB8(J) ,J=I, 3) 

IA  =  TAB8(  1 )  *TA88(2)*TAB8(  3M-3.0 
REAC  (5,3112)(TAB8(J),J=4,IA) 


RE  AC  (5,3112) (TA89{ J),J=i, 3) 

I A  =  FA89 (  1 )  *TA89( 2) *TAB9( 3) +3.0 
READ  (5,3L12)(rA89(J),J=4, 1A) 

2002  IF  (R. Nt. LANGUR. AND. LANGMR.NE.O. )  GO  TO  2005 
REAR  (5,3100)  HEAD 
IF  ( FEAC ( 1 ) . EC. I  BLANK )  GO  TO  2004 
RE  AC  (5, 1000) ( r A { I ) , 1  =  1,6) 

READ  (5, 1000) (T( I ), 1*1,6) 

T  A ( l ) =-23. 

RE AC  (5,1000)  ALWC,R,UINF,UPINF,PA,WA 

REAC  (5,1000)  X.MAX,YMAX,YDROP»DELYOP»XL»XHLi 

RE AC  (5,1000)  YSL,YSY,VXWIN,VYWIN, DIRECT, RANDOM 

REAC  ( 5 , ICOO )  XHL, YHL , Al AMDA , RE S I L , RHOW , t ANGMR 

REAC  (5,1000)  GRAVTY,CELX1,DEI.X2,XREF1,XREF2 

REAC  (5,1000)  YTRAP1,YTRAP2,XTRAP1,XTRAP2,TABLE 

SY0R0P=YDR0P 

I BEV=  1 

NK=0 

KL  =0 

IF  ( L ANGMR. NE. 0. )  R=L ANGMR 
DO  2CC7  1=2,300,2 
2007  rA811(I)=0. 

I MAX=0 

2C05  REl.\'F  =  2.0*R*PA*UINF/WA 

IF  (FABLE  .EG. 0.0)  GO  TO  6055 
T  A  8  1  (  1 )  =  l , 
r  AB 1 ( 2 ) =2  • 

T  AB 1 ( 3 )  =  100  • 

TA8l(4)=0, 

T  A8 l ( 5 ) =0  » 

NK  =  5 

00  701  1=6, 100,2 

7C l  IF  (I A86( I ) .G6.XML)  GO  TO  702 

702  NK6  =  I - 1 

703  NK=NK+1 
NK6=,MK6+1 

r  A8 l { NK ) =TAB6(NK6)-XHL 

NK=MK+1 

NK6=NK6+l 

IF  (TABl(NK-l).EC.O.  )  GO  TO  705 

TAlU  (NK)=TA81(NK-2)  +  SQRT( ( TAB6 ( NK6- 1 ) -TAB6 ( NK6-3 ) ) *  *  2+ ( T  AB6 ( NK6 ) - 
XTAUMNK6-2) ) *  *2 ) 

IF  ( r  ABM  NK6-1 ) .LT .XMAX-.COl )  GO  TO  703 
GO  FC  706 
7C5  TAB l ( NK ) =0. 

GO  TO  703 
7C6  TA821 l  )  =  i. 

NK=(NK-3)/2 
T  AB l ( 3 ) =NK 
FAB2(2)=2. 

T  A  8  2 ( 3 ) =100* 

T  A  8  2  ( 4  )  =  0  • 
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r A82 ( 5 ) =0, 

NK  =  5 

CO  7 10  1=6,100,2 
IF  (TAB7m.GE.XhU  GO  TO  712 
NK7=!-1 
NK=NK+l 
NK7=NK7U 

TAB2 ( NK ) =TAB7 (NK7 l-XHL 
NK=NK*i 
NK7=NK7  + 1 

IF  (TA32(NK-l).EQ.O.)  GO  TO  715 

TAB2(NK)=TAB2(NK-2)+  <0R  T  f ( TAB7 ( NK7- 1 } -TAB7 ( NK7-3 ) ) * *2+ ( TAB7 ( NK7 ) 
XT  A 157  { NK7-2 )  )  **2 ) 

IF  (TA87(NK7-i) .LT.XMAX-.OOOl )  GO  TO  713 
GO  TO  716 

715  T AB2 ( NK ) =0. 

GO  TO  713 

716  T AB9 { 1 }  =  1 • 

T A89 { 2  J  =2 . 

TAB 9(31=50 
NK=(NK-3)/2 
T  AB2 ( 3 ) =NK 

T A 39(41=0. 

T  AB9 ( 5 1 =0  • 

NK  =  5 

CO  720  1=7,50,2 

720  IF  (TAB8( Il.GT.O.)  GO  TO  722 

722  NKfl  =  1-4 
XHLCB=TAB8( 1-3) 

723  NK=NK+i 
NK  8  =NK8  +  l 

TAH9(NK)=TA88(NK8  l-XHL 
NK=NK+  1 
\:<0  =  NK8  +  1 

IF  (NK.6Q.71  GO  TO  725 

T  A39 { NK 1 =TA89 ( NK-2 1 +SQR  T ( { TABS ( NK 8- 1 1  - TAB8 ( NK8-3 1 1  *  * 2 M  T AB8 ( NK8 1 
x  T  A33 ( NK8-2 1 1**2) 

IF  ( TA38 ( NK8- 1 ) . IT. X VAX-. 0001)  GO  TO  723 
GO  TO  726 

725  T  A3'9  (  NK  )  =0. 

GO  TC  723 

726  NK  = ( NK- 31/2 
T  A  3  9  ( 3 ) =NK 

60  5  5  I  A  =  r  A 8 6 ( l 1  * T AB6 ( 2 1 « T A B6 ( 3 1 
WR I T  t  (6,670) 

!A=TABU 1)*TABI(2)*TAB1( 31+3.0 
WRITE  (6,680)  ( TAB1 ( I  1 , I =6, I  A) 

WRITE (6, 605) 

WRITE  (6,671) 

I  A  =  T A82 ( l ) *TAB2( 2)*TA92( 3) *3.0 
WRITE(6,680)  ( T  AB2( l ) , I =6, I A ) 

WR I TE ( 6 ,605  ) 


710 

712 

713 
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wr i re ( 6 , 67  3 ) 

I  A=  TAB4 ( 1 )*TA34{ 2)*TA8M 3) *3.0 
WRITE  (6,680)  ( TARA (  I )  ,  l =6 , I  A ) 

WRITE  (6,605) 

WRITE  (6,674) 

I  A=  T A85 C 1 ) *T AS5 { 2 ) *  TA85  C 3) +3.0 
WRITE  (6,680)  ( TA85(  I  ) , 1=6, I A ) 

WRITE  (6,605) 

WRITE  (6,675) 

I  A  =  TAB6  {  1 ) *TAB6( 2)*TAB6( 3) *3.0 
WRITE  (6,680)  { TA86( I ) , 1=6, I  A ) 

WRITE  (6,605) 

WRITE  (6,676) 

1  A  =  TA87 l I)*TAB7(2)*TA87{3)+3.0 
WRiTE  (6,680)  ( TA87 ( I ) , I =6 , I A ) 

WRITE  (6,605) 

WRITE  (6,677) 

I  A  =  T A 8 8 { I ) *T A88 ( 2 )*TAfi8 ( 3 ) »3,0 
WRITE  (6,680)  (TA88U)»  1=6,  IA) 

WRITE  (6,605) 

WRITE  (6,678) 

I  A  =  T A 8 9 { I ) *TAB9( 2)*TAB9{ 3) *3.0 
WRITE  (6,680)  { TA89 ( I ) , I =6 , I  A ) 

WRITE  (6,605) 

T  ABLE =0.0 
Y0RGP=SY0R0P 

AK=2.0»Rh0W»R*«2*UINF/(9.0»WA) 

WRITE  (6,3iOL)H£AO 
WRITE  (6,3102) 

WRITE  (6,3106)ALWC,R,UINF,UPINF,PA,WA 
WRITE  (6,3108) 

WRITE  (6,31065  XMAX, VMA X , YDR0P , DEL YOP , XL , XHl i 
WRITE  ( 6 , 3  L  0  4 ) 

WRITE  (6,3106)  YSL,YSM,VXWIN»VYWIN, DIRECT, RANDOM 
WRITE  (6,3105) 

WRITE  (6,3106)  XHL  ,  YhL  ,  ALA.MOA  ,  RES  I  L  ,RHOW  , LANGNR 
WRITE  (6,3114) 

WRITE  (6,3106)  GRAVTY,DELX 1 , DELX2 , XREF i , XREF2 
WRITE  (6,3115) 

WRITE  (6,3106)  YTRAPl,YTRAP2fXTRADl,XTRAP2,TA8L£ 

WRITE  (6,3165) 

,,’ITG  (6,3106)  RE  INF  ,  AK 
WRITC  (6,3107) 

00  2010  1=1,500 
20  10  T  A8  10  (  I )  =0 . 

T  A  S  1 0 (  1  )  =  1. 

TAB10(2)=2. 

TAR1CI 3)=500. 

TA  *3  10(41=0. 

TA310(5)=0. 

xpy 1=0. 

XXT  =  0.0 
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CX  *  0.0 
IP  =  0 
I  PA  =  0 
JP  =  0 
YP=0.0 
YPM1  =  C.o 

watrat  =  o.o 

AVS  =  0.0 
UNDER= 1 . 0 
VXW=VXWIN 
V YW=VYW IN 
STRCIR=CIR£CT 
FN=0.0 
PM.M=0.0 
GO  ro  104  ' 

LOG  YM=YKNi 
YN=ABS { YH) 

XNN* XNN L 

HUNX\N= 1 00 . *  XNN 1 
IF(T(5)  -XNN) 3400, 3400,3404 
34CC  I  F ( T ( 6 )  -  YM) 3401,3401, 3404 
3401  I F  {  T {  I )  -XNN)3404, 3402, 3402 
34C2  I  F ( T ( 2 )  -  YM  )  3404, 3403, 3403 
3403  JP  =  JP  +  L 

JP3  =  IPA  JP 

IF  ( [PB.LT.3099)  GO  TO  6064 
WRITE (6, 662) 

I Pfi= I  PA 
•)P=0 

6064  XPLTI IPQJs  XNN 
YPLT(  I  Pt3 )  =  YM 

3 4 C 4  10  *  ID  M 

\ 

IF  (tC.LT.549)  GO  TO  6065 
WRITE  16,661) 

io=tc-i 

6065  XPLCTUO)  =  XNN 
YPLCTUO)  =  YM 
VOXW=VOXWNi 
VOYW=VOYWNl 

GO  TO  2205 

IC2  IF  (RANCOM)  305,305,310 
3C5  YDRCP=YQROP+OELYDP 
3C6  IF  (RANDOM. GT. 0.0}  GO  TO  310 
VX'W  =  VXW I N 
V YW  =  V YWIN 
GO  TO  315 

310  R  £ AC  (5,1000)  YCROP , XL  >  VXW , VYW 
3CUNCE=0.0 

IF  (YOROP+XLfVXW+VYW.EQ.O.O)  GO  TO  2003 
IF  (YCROP. GT.YMAX)  STOP 
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315  I F  {  J  P }  3405,3406,3405 
3^05  ! P  =  IP  +  l 
I  PC l  IP)  =  JP 
I  PA  =  IPA  +  JP 
JP  =  0 

3406  T A { 7 }  =  ID 

DIRECT=STRD!fi 

CHI  ^PLOTStTA,XPLOT,YPLOn 
104  VfJXw-VXW/UINF 
V 0YW=V YW/U INF 

WRITE  {6,603}  XL  *  YDRGP,  VOXW,  VOYVf 

Yy=YCROP 

YN=ABS(YH) 

WRITE  (6,3109) 
hUMXNN=IOO,*XL 
XNN=XL 
XNN 1= XNN 
XPLCHl)  -  XNN 
YPLGTCU  a  YH 

iF(TfS)  -X\'N}3407r34C7,34U 
340?  I F l T { 6 }  -  YM I  3408  *  3408, 3  4 1  £ 

3408  i?(T{i)  -XNN ) 341 l , 34C9, 3409 

3409  I  F  {  T  {  2  )  -  YM  j  34 i L , 3 4 i 0 , 34 1 0 
34 1C  JP  =  JP  *  i 

l  P8  =  IPA  4-  JP 
XPLTUPB)  =  XNN 
YPLT (  I  P8  }  =  YM 
3411  ID  =  2 

XPLCniO)  a  XNN 
YPLOT { ID)  =  YM 
IF  rXNN.GT.XU  GO  TO  2205 
103  \YR*TA83(2)-l. 

NXC*TA83<3)-1. 

CO  210  I=i,NYR 

Y( I ) =FLOAT ( I  NT { T  Ad  3 ( J  +4 } *  10000. 0  +  . 5 ) )/ 10000.0 
T A 8 3 ( I+4)=Y(  I) 

2LC  CONTINUE 
J  =  4 

DO  220  I  =  1 » NXC 
J  =  J  +NYR  +  l 

X  (  t  ) -FLOAT  {  I  NT  (  T  A03  (  J  )  »  l.COOO .  0+  .  5  )  )  /  10000 . 0 

T  A  8  3  (  J  )  =  X  { I  } 

22C  CONTINUE 

2205  IF  (YN.GT.YKAX)  GO  TO  2003 
DELX«CELXl 

IF  (XNN.GT.XREFl .AN0.XNN.LT.XREF2 )  0ELX=DELX2 
221  IF  (CIRECT)  225,225,260 

225  IF  (VCYW)226, 226,235 

226  IF  (YM.LE.O.)  GO  TO  231 
232  CO  2 28  1=1, NYR 

IF  (Y(I ) .GE.YN)  GO  TO  229 
228  CONTINUE 


229  Y  M  N 1  =  Y  (  1-1) 

2  3C  Y J  =  Y ( 1-1) 

IFIYN.LE.-YI  1)  )  YMN1=-YMN1 

IF  {  A6SIYNN15.EO.YI2)  .ANO.VDYW.GT.O.O)  YM.NI  =  AQS  I  YKNl ) 
YJ1  =  Y {  I  ) 

GO  TG  250 

235  IF  IYP.LT*G.)  GO  TO  232 
231  CO  236  I  - 1 f  NYR 

IF  miKGT.YN)  GO  TQ  237 

236  CONTINUE 

237  YMN I =Y { I  ) 

GO  TO  230 

250  00  255  1=1, NXC 

IF  { X (  I )  .GT.XNN)  GO  TO  256 

255  CONTINUE 

256  XN1  =  XU) 

X  N  =  X  (  I  —  L  ) 

IF  (CIRECT.GT.O.  )  XNN1  =  ( HUNXNN+OELX ) / 100. 

GO  TO  13 

26C  CO  265  1=1, NYR 

IF  (YU  )  oGE.YN.)  GO  TO  266 

265  CONTINUE 

266  IF  (YN.EO.Y(D)  1=2 
YJl=Y(  I  ) 

YJ=YII~!) 

IF  (VCXW)  270,271,250 
271  IF  (VOX)  270,270,250 
27C  CO  275  1=1, NXC 

IF  (XII ) .GE.XNN)  GO  TO  276 

275  CONTINUE 

276  XN 1  =  X (  I  ) 

XN  =  X I  I  - 1  ) 

X  \  N  l  =  I HUNXNN-OELX ) /  1 00 . 

13  CALL  rLUXXSI TAB6,XNl,0,0,YB2) 

CALL  TLUXXSI TAB7,XN1,0,0,YB3) 

YB2NI  =Y82 

Y  P  3  N  J.  =Y83 
IMO 1  =-l 
INC2  =  -2 
I  NO 3  =  -3 
I N04=-4 

IF  { XNN1-XHL I  )  11,12,12 
12  I F I YN-YSL  )  11,1301,1301 
1301  I  F  ( YN-YSN* )  1 A ,  14,  11 

14  CALL  TLUXXSITA34,  YJ,0»0,XB2) 

CALL  TLUXXSI TA85,  YJ,C,0,X83) 

XR2 J=XB2 

XB3J=X83 

CALL  TLUXXS(TAC4,YJ1,0,0,XB2) 

CALL  TLUXXSI TA85,YJl ,0,0, XB3) 

XR2J1=X82 
X83J l =XB3 
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CALL  TLUXXS! TAB6 , XN , C , 0 , Y3 2 ) 
CALL  TLUXXS!  7 AB7  ,  X,I ,  C , 0  ,  YB 3  ) 
YB2N=YB2 
YB3N=YB3 

CALL  TLUXXS(TAB6,XN1,0,0,YB2) 
CALL  TLUXXS!TAB7,XN1,0,0,YB3) 
CALL  TLUXXS  (TA88,XN,0,0,YBM 
YB1N=YBI 

CALL  TLUXXS  (TAB8,YNl,0,0fYBi) 

YB1N1=YBI 

YB2N'1  =  YB2 

YB3N 1 =Y83 

IND1  =  -1.0 

I ND2=  - i . 0 

I  NC3  =  -1.0 

INDA  =  -1.0 

IXB2J=ICCC0.0*(XB2J+. 00001) 

I X83 J = 10000.0* I XB3J+. 00001) 
IXB2J1=10000.0*(XB2J1+0. 00001) 
IXB3JI=10C00.0* ( X 83 J 1+0.00001 ) 
IYB2N  =  10000. 0MYB2N  +0.00001) 
IYG3N  =10000. 0MYB3N  +0.00001) 
IY82N 1=10000, C+(YB2Nl+0.C000i) 
IYB3'N-1  =  10000. 0*(Y33NJ  1+0. OCOOl) 
IYBIN  = 10CC0 . 0* ( YB IN  +0.CC001) 
I Y8 1M1  =  10000. 0*(YB IN  1+0, 00001) 
IX, N  =10000. 0*(XN  +0.GC001) 

IXN1  =  1 OCOO  .  GM  XN  1  +0.00001) 

IYJ  =1CC00.0*(YJ  +0.00001) 

IYJ1  =ICC00.0*IYJ1  +0.0C001) 

IF  (  l  XN-  I  XB2 J )  15,16,16 

15  X  L  L  =  X  N 
YLL  =  YJ 

20  IF  { I XN-  I X32 J  l )  17,13,18 

16  IF  ( IXN.GE. I XB3 J ) GO  10  15 
IF  ( IYJ. GT. IYBIN)  GO  TO  19 
XLL=XB3J 

YLL  =YB IN 
INC! =2 
GC  TC  20 

19  XLL=XB3 J 
YLL  =  Y  0  3  N 
I  NO  1 =  1.0 
GO  TO  20 

17  X  U  L  =  X  N 
YUL  =YJ  1 

22  IF  UXN1-IXB2J)  2  3  »  2  A  ,  2  A 

13  IF  (  IXN.GT. IXB3J l )  GO  TO  17 
IF  (  IYJ. GT. IYBIN)  GO  TO  21 
XUL=XB3J 1 
YUL=Y81N 
IN02  =  2 
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GO  TO  22 
21  XUL=X83JI 
YUL=YB2N 
I N02=  1.0 
GO  TO  22 

23  XLR  =  XN 1 
YLR=YJ 

281  IF  { IXN1-IXB2J1 )  25,26,26 

24  IF  ( IXN1.GT. IXB3J)  GO  TO  23 
IF  (  IYJ.GT.  IYBINDGO  TO  27 
XLR=XB2 J 

YLR=YB1N1 
INC3=2 
GO  TG  281 

27  XIR=X82J 
YLR=Y83N1 
INC 3=  1*0 
GO  TO  281 

25  XWR=XNl 
YWR=YJ1 
GO  TO  29 

26  IF  { IXN1.GT • IXB3J1)  GO  TO  25 
IF  ( IYJ.GT. IY81N1I  GO  TG  28 
XKR=XB2Jl 

YWR  =Y8 INI 
INC4=.2 
GO  TO  29 

28  X*R=X82J1 

Y  V<  R  =  Y  8  2  N  l 
INC 4=  1.0 
GC  TO  29 

11  YLl =Y J 

Y  U  L  =  Y  J  1 
YLR=YJ 
YWRsYJl 
XLL=XN 
XUL  =XN 
X  L  R  =  X  “  1 
XWR=XN1 

29  C£L YL=YUL-Yll 
IF(CELYl)  30,30,31 

30  VXL  =0 

32  CELYR= YWR-YLR 

IF(OtLYR)  33,33,34 

31  IF(  INCl  )  43  1  ,  43  1  ,  432 

432  PHlll  =  1.0 

IF  IINC1.EQ.2)  PH  I It  =0 . 

GO  TO  433 

431  CALL  TLUXXS ( TAB3 , XLL  »  YLL , 0 ,  PHILU 

433  IF ( INC2)  434,  434,  435 
435  PHIUL  =  1.0 

IF  (INC2.EQ.2)  PH IUL  =0. 


I. 
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GO  re  46 

434  CALL  TLUXXSITA83,  XUL,  YUL ,  0,  PHIUL) 
46  IF  ( ALAPDA.6G.0.0)  GO  TO  436 
DEL YL= ( I YLL^YUL )  /  2. ) *OELYL 

436  VXL=(PHIUL-PHILL)/OELYL 
GO  TO  32 

33  VXR=0 

35  CONTINUE 

3221  OELXL  =  XLR  -  XLL 
I F { CELXL  )  36,36,37 

34  I F ( IND3)  437,  437,  438 
438  PHILR  =1.0 

IF  (INC3.EQ.2)  PHILR=0. 

GO  TC  439 

437  CALL  TLUXXSI TAB3 , XLR , YLR ,0 ,  PHILR) 

4  39  IFUNC4J  440,  440,  441 

441  PHIWR  =1.0 

IF  (1NC4.EG.2)  PHIWR=0. 

GO  TC  3010 

44C  CALL  TLUXXSI TAB3 , XWR , YWR , 0 , PH  I WR ) 

30 1C  IF  { ALAMDA.EG.O.O)  GO  TO  442 
CEL YR= ( (YWR+YLRJ/2. )*C£LYR 

442  VXR  =  (PHIWR-PHILR)/0GLYR 
GO  TO  35 

36  VYL  =  0 . 0 

38  OELXU=XWR~XUL 

IF(CELXU)  39,39,40 

37  I F ( I NC3  )  443,  443,  444 

444  PHILR  =1.0 

IF  (INC3.EQ.2)  PH  I LR  =  0. 

GO  TO  445 

443  CALL  TLUXXSI TAB3, XLR, YLR, 0, PHILR) 

445  I F (  I  NO  l )  446,  446,  447 

447  PHILL  =1.0 

IF  ( I N  C 1 » £  Q ♦ 2  )  PH  ILL =0* 

GO  TO  448 

446  CALL  TLUXXSI TA83, XLL, YLL,0, PHILL) 

448  IF  I ALAMOA.EG.O. )  GO  TO  4481 
DELXL=Y J*DELXL 

4481  VYL=  (PHILR-PHILL 3/OELXL 
GG  TC  38 

4 C  IF!  INC4)  449,  44°,  450 
45C  PHIWR=  1.0 

IF  (INC4.EQ.2l  PH  I WR  =  0. 

GO  TC  451 

445  CALL  TLUXXSI TAB3, XWR , YWR, 0,PHIWR) 

451  IF! I N  0  2  #  452,  452,  453 

453  PH  I Ul.  =  1.0 

IF  UNC2.EQ.2)  PHIUL=P. 

GO  TC  454 

452  CALL  TLUXXSI TAB3, XUL, YUL, 0, PHIUL) 

454  IF  ( ALAMOA.EG.O.  )  GO  TO  4541 
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DELXU=CELXU*YJl 
4541  VYW  =  { PHI WR-PHIUl i/OELXU 
WFS  =0 . 0 

co.nter=o.o 

GC  TC  4101 
39  VYW=0 .0 
WFS =0.0 
C0NTER=0.0 
4101  CIVX=0.0 
0  I VY  =  0  «0 

IF  (CIRECT)  320,320,340 
320  IF  (VYW.NE.O.)  C1VY=1. 

IF  (VYL.NE.O.)  0  I VY  =  C  l  VY  +  1 * 

IF  { YXL  )  360,361,360 

361  VX=VXR 

GG  TG  364 

36C  IF  { VXR }  363,362,363 

362  VX=VXL 

GO  TO  364 

36  3  V  X  =  V  X  R  *  ( (XNN-XLL)/( XLR-XLL) ) + VXL  * ( { XLR-XNN ) / ( XLR-XLL ) ) 

364  YF'EAN  =  AB${  ( Y MiM  1  +  Y M }  / 2 . 0 ) 

VY=(VYLMYP6AN-YJi ) +VYW* ( Y J-YMEAN ) ) / ( YJ 1-YJ ) 

GO  TC  3011 

340  IF(VXL)  42,43,4 2 

42  DIVX=1.0 

43  IF(VXR)  44,45,44 

44  DIVX-OIVX  +  UO 

45  X.v  FAN  =  (  XNN+XNN1  )/2. 

VX=(VXL*(XN1-XY6AN)+VXR*( XME AN-XN ) ) / ( XNl-XN  J 
IF  ( VYW.NE.O.O.ANO.VYl.NE.O.O)  GO  TO  365 

IF  (VYW.NE.O.)  GO  TO  7001 

VY=-VYL 

GC  TC  3011 

7C0 1  IF  (YJ.EQ.O.  )  GO  TO  365 
VY=-VYW 
GO  TO  3011 

365  VY  =  -VYW* ( ( YN-YLL ) / { YUL-YLl ) )-VYl*((YUl-YN)/(YUL-Ytt) ) 

3011  VOX  =  VX/UPINF 

VOY= VY/UP I NF 
IFIYN'.LT  .0»  )  Vi)Y  =  -VDY 
F  IRST  =  0. 

0WB  =  0 .0 
VDXWM=VCXW 
VCY  ,Viv  =  VOYW 

325  RC\=KEINF*SQRT{ { VOX- VOX WM ) « * 2  +  ( VDY- VOY WM ) » *2 ) 

CK=l.+. 19 7*REN**. 63 +.00026 *R£N*» 1.38 
IF  (DIRECT)  336,336,337 

3  36  XX5X  =  VDYW**2  +  (YM,Ni-YM) * ( VOY-VDYWM) *CK/ { 6 . 0 * AK ) -5 . 3623/ (RHOW«UINF»» 
X  2 ) * ( RFCW-PA ) * ( YMN 1-YM ) *GR  4 VTY 
GO  TO  338 

337  XX5X=VDXW«*2+( XNN1-XNN) * ( VDX-VDXWM) *CK/(6.0»AK) 

339  IF  ( XX5X.LT .0. )  GO  TO  370 
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If  (WFS.GT.O.O)  XX5X=0.0 
IF  (DIRECT.GT.O. )  GO  TO  339 
VDYWNL=SCRT{XX5X ) 

IF  (VCYW.LT.O.)  VCYWN i=-VDYWNl 
0rw=(YMNl-YH)/(6.C»( VOYWN1+VOYW) ) 

VOXWNl=VCXW+CTW» ( VOX-VDXWM } *CK/AK 
GO  TO  5010 

37C  IF  (CIRECT.GT.O.O)  GO  TO  348 
IF  (CCNTER.GE.4.0)  GO  TO  348 
C0NTER=C0NTER*1.0 
Ym  =  (  YM*YMNi)/2.0 
GC  VO  4101 

348  IF  (CWB.GT.O.O)  GO  TO  345 
OWBsl.O 

IF  (DIRECT. GT. 0.0)  GO  TO  346 
V0YWiv  =  VDYW/2.0 
GO  TO  325 

346  V0XWM=V0XW/2.0 
GO  TO  325 

345  IF  (DIRECT. GT. 0.0)  GO  TO  347 
C I RECT  =  -C  IRECT 
GO  TO  2205 

347  XNNl=XNN-(VDXW*VCXW»6.0*AK)/(CK*{VDX-V0XW/2.0) ) 

conter=conter+i.o 

W  F  S  =  1 . 0 

IF  (CCNTCR.LE. 10.0)  GO  TO  4101 

VOYWf/  =  VOYW 

XKT=26.0 

XX5X=0.0 

GC  TO  339 

339  V0XWN1=SQRT{XX5X) 

IF  (VCXW.LE.O.)  VCXWNls-VDXWNi 

IF  (VCXW.EQ.O.O.ANO.VDX.GT.O.O)  VOXWN 1 -AOS { VOXWN 1 ) 

DTW=  {  XMN  l-X.NN  )  /  (  6 .0  *  {  VOXWN  1  +  VDXW  )  ) 

VDYWN1 =VOYW  +  C  TW* { VOY-VOYWM ) *CK/AK-32.  1 74*DTW * { RHOW-P A ) / ( RHOW*U I NF * 
X  »  2 ) »GRAVT Y 

30 1 C  IF  ( VCXWN! .GT.VCXW)  GO  TO  4200 
IF  ( VDX.GE.VDXWNI )  VCXWN1=VDX 
GO  TO  5020 

4 20C  IF  { VCX . L E . VCXWN  1 )  VCXWN1  =  VDX 
502C  IF  ( VCYWNi.GE.VUYW)  GO  TO  4100 

WDY  =  V0Y-32.17*»*(  RhQW-PA  )  »  {  DTW  )  /  (  RH0W»UINF*»2  )  *GR  A VTY 
IF  ( VCYWN l »GE . WDY  )  GO  TO  399 
VDYWN 1 =VOYW 

IF  (WCY.LT.VCYW)  VOYWNl=WOY 
GO  TC  399 

4  ICC  WOY  =  VDY-32. 1  74*{RHOW-PA)» ( DTW ) / ( RHOW*U I NF **2 ) »GR AVTY 
IF  (WCY.GT.VCYWN1 )  GO  TO  399 
VDYWN l =VDYW 

IF  (WOY.GT.VCYW)  V0YWN1=WDY 
399  V0XWK1  =  { VCXWN 1  ♦  VCXWI/2.0 
V0YWM1  =  (VDYWN  1  *  VCYWJ/2.0 
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I  F  C  X  K  T  -  25.0)4999,4*399,3028 

4999  I F { A fi S  C VDXW )  -  .01  )3041,3040,3040 

304C  IFMBSfVCXWM  -  V0XKM1)  -  ABS  (  .001  *  VDXWM )) 304 1 » 304 1 , 3026 
3041  I  F  (  A6S  ( VDYW )  -  .01)3028,3024,3024 

3024  I  F  (  ABS  ( VOYWM  -  VCYWM1)  -  ABSC.OOl  »  VOYWM )) 3028 , 3028 , 3025 

3025  I  F  (  ABS  ( VOYWM  -  VOYW'Ml)  -  .CC01  )  3028 , 3028 , 3027 

3026  l F ( ABS (VDXWM  -  V0XWM1)  -  .0001  ) 304  1 , 304 1 , 302 7 

3027  OVXV=VCXWM1-VOXWM 
CVYM=VOYWMl-yDYWM 

IF  (FIRST. L£. 0.  )  GO  TO  750 

VOYwM 1 =VCYWM-DVYMM  {  VDYWMX- VDY >. V  )  /  {  DVYMX-DVYM )  ) 

V cxwyi* VOX WM-DVXM *( VCXWMX- VDXWM) / ( DVXMX-OVXM ) 

75C  VDXWMX=VDXWM 
VDY'rtMX=VDYWM 
VCXWM=VDXWM1 
VOYWM=VOYWM1 
FIRST*FIRST+1. 

CVXYX=CVXM 
CVYMXaOVYK 
XKT  =  XKT  ♦  1.0 
IF(XKT-25. ) 325 , 325, 5000 

5000  WRITE  (6,3111 ) VOYWM , VDXWM 
voxwy  =  vOXW 

VDYh'M  =  VOYW 
GG  TC  325 

3028  IF  (CIRECT)  500,500,502 

5CC  XNN  1  =  XN\'  +  6 . 0*DTW«  {  VDXKN  l+VDXW ) 

XKT  =XK  T+  l  .0 

IF  ( ABS ( XNN 1-XNN ) .L6.0ELX/55.0)  GO  TO  52 
46C  IF  (XKT. GT. 26.0)  GO  TO  52 
DIRECT  =-DlR6CT 
GO  TC  2205 

502  Y.v\  l  =  YM*-6.0«CTW»  (  VOYWNi+VOYW) 

XKT  =  XKT  +  1 .0 

IF  (ABS(YM.M1-YM).LE.(YJI-YJ)*1.05)  GO  TO  52 
GO  TC  460 

52  WRITE  (6,3110)  XNN 1 ,  YMN  1 » DTw ,  CK  ,  REN  ,  VOX  ,  VDXa'i  l »  VD  YWN  l  ,  VDY ,  XKT  , 
XDIRECT 

XKT  =  0.0 

3020  IF  (XNNI-XMAX)54,53, 53 

53  IF  (ABSCYWN1  )-YSM)  459,459,505 
459  IF  (XPMl.EQ.O.O)  GO  TO  462 

IF  (FNcNE.0.0)  GO  TO  471 
SAVE*YCROP 
F  N  =  1 . 0 
GO  TO  461 
471  FN=FN+1.0 

461  YCR0P=Y0RCP-CELYCP/2.0*»FN 
GO  TC  306 

462  IF  (FN.EQ.O.O)  GO  TO  102 
FN=FN+ 1.0 

YDR0P=YDR0P+0ELYCP/2.0»»FN 
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GO  TC  306 

505  IF  (RANDOM)  5070,5070,102 
507C  IF  (LANGMR.LE.O. }  GO  TO  2003 
ND=  5  +  2*  t NK- 1 ) 

WRITE  (6,630)  R, (TA310(I), 1=6, NO) 

AVS'VAX  =  TA610{N0-I  ) 

AVS=FLCAT { INT{TA810{6}*4.}-~i)/4. 

GO  TC( 5040, 5041, 5042,5043, 5044,5045,5046) , IBEV 

5040  FACTOR3 • 3 

5048  APERAG-.71 
GO  TC  5055 

5041  FACTOR3. 2 

5049  APEQAO3 . 52 
GO  TC  5055 

5042  FACTORS. I 
5C5C  APERAC-. 31 

GC  TO  5055 
5C43  F  AC  TGR= . 05 

5051  APERAC=2.22 
GO  TC  5055 

5044  FACTOR3. 05 

5052  APERAG= 1.74 
GO  TC  5055 

5045  FACTCR=.i 

5053  A p £ R AO  = 1 . 3 7 
GO  TC  5055 

5046  FACTOR3. 2 

5054  APERAC=0.0 

5055  R  =  APERAO*LA,NGMR 
BEV  =  I BEV 

/<R  I  TE  (6,640)  BEV, FACTOR,  APERAO 
IBE7=I8EV+l 

:>G 56  C  ALL  T LUX XS  ( TAB  10 , AVS , Q ,0 , WATRAT ) 

W  RI  T  E  (6,640)  WATRAT 

IF  ( AVS.GT.AVSMAX)  W  A  T  R  A  T  =  0 . 

WAWRT=FACT0R*WATRAT 
00  5072  1=1,800,2 

5 C 7 2  IF  (TABiim.EQ.AVS)  GO  TO  5073 
5073  CUMWAT  =  TABl  1  (  1  +  1) 

WRITE  (6,640)  CU^WAT  ,WAWRT , AVS 
CUMW  AT  =CUMW AT  +WAWRT 
TABIK  1  +  1  )=CUMWAT 
AVS=AVS+0.25 

IF  (V.ATRAT  .GT.O.C)  GO  TO  5056 
I  E\Cr*2*NK  +  5 
00  5C60  1=6,  I  ENO 
506C  T  A3  10  (I)=0. 

N  K  =  0 

IF  (R.LE.O.)  L ANGMR  =0 . 

IF  (I.GT.1MAX)  IMAX= I 

WRITE  (6,635)  ( T  AB  1  l  ( J  )  , J=  l ,  IMAX  ) 

GO  TC  6041 
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54  IF  ( ABS ( YMN 1 J-YMAX )  55 »  53  »  5 3 

55  CALL  TLUXXSITAB6,XNN1,0,0,Y82) 

WFSMB=0. 

XNI=XNN1 

IF  (ABS  (YMND.Lr.YB2)  GO  TO  L575 
CALL  f LUXXS ( TAB7 , XNN  1 , 0 , 0, Y83 ) 

IF  ( ABS ( YMN 1 ) -YB  3 )  57,100,  100 
1575  CALL  TLUXXS(TA03»XNN1»0,0»YB1) 

IF  {  A  8  S  ( YMN  D.GE.YB1)  GO  TO  100 
1577  CALL  TLUXXS  ( T ABS , XN  1 , 0 , 0 , YB D 
YB1N1=YB1 

CALL  TLUXXS(TAB8,XNN,0,0,YB1) 

YB1N=YB 1 

IF  (  XNN. GE. XhL  1 }  GO  TO  472 

Y2 1  * ( YM* ( XNN1-XHL 1 ) +  YMN 1  * ( XHL i-XNN )  ) / { XNNl-XNN ) 
AMQ=Y81N1/(XN1-XHL1) 

AMCT= ( YMN1-YM ) / ( XNNl-XNN) 

IF  IYZ1.GT .0.0)  GO  TO  473 

CIFX=(0.0-YM-AM0T*(XHL1-XNN)  )  / t  AMB  +  AMDT ) 

XP=XHL1+DIFX 

YP=-AMB*OIFX 

AMC  T  =  -AMOT 

GO  TC  474 

4  73  G I F  X  = ( YM  f  A  M  0  T  * { X  h  L 1 - XNN )  )/(  \M3-AM0T) 

XP=XhLl+OIFX 
YP=AM8*0IFX 
GO  TO  474 

4  72  A  M  B  = ( Y  B 1 N  l  - Y  8 1 N ) / ( X  N  1  - X  N N ) 

AMOT =( ABS ( YMNl )-YN)/( XNNl-XNN) 

0  I F X  = ( YN-YB IN ) / ( AMQ-AMDT ) 

XP  =  X.\N  +  OIFX 
YP=YB1N+ AMB'DIFX 
IF  (YMN1.LT.C.0)  YP=-YP 
4  74  CON T  I.  NUb 

IF  { ABS ( XP-XN i  ).LE. 0.0001)  GO  TO  3501 
X  N 1  =  X  P 

W  F  S  M  8  =  W  F  S  M  8  ♦  l  •  0 
i F  (WFSVB.LE.25.0)  GO  TO  1577 
GO  TC  3501 

57  CALL  TLUXXS  { T A86 , XN 1 , 0 , 0 , YB 2 ) 

Y62Nl=YB2 

CALL  TLUXXS  ( T  A  B  6 , X  N  N , 0 , 0 , Y  B  2 ) 

Y  B  2  N  =  Y 13  2 

IFIYB2N-YMAX)  157,158,153 

158  IF  (XHL.LT.XN1)  GO  TO  1585 
XNI-XNNI+CELX2/100. 

GO  TC  57 

1585  YZ= ( ABS( YMNi ) M XhL -XNN) +YN* { XNNl-XhL ) )/( XNNl-XNN) 
IF  (ABSIYZ ) -YHL »  157,3500,160 

159  AV8=(YB2N1-YFL)/(XN1-XHL) 

GC  TC  161 

16C  CALL  TLUXXS  ( T AB 7 , XN 1 , 0 , 0 , YB 3  ) 


Y83M=YB3 

AMR = { YB3N 1 -Ybl ) / ( XN 1— XHL ) 

161  Y0=YFL-AP8» ( XHL-XNN) 

GC  TO  60 

157  I F { YN-YI32N )  58,  59,  59 

58  A N* p  =  f  YB2N1-Y32N)/I XNl-XNN) 

CALL  TLUXXSt TAB6,XNN,0,0,YB2) 

Y0=Y82 

GC  TC  60 

59  CALL  TLUXXS  { TAB7 , XN l , 0, 0, Y83 ) 

Y83N i=YB3 

CALL  TLUXXS(TAB7,XNN,0,0,Y63) 

Y0=YB3 

APB={YB3\i-Y0)/ I XNl-XNN) 

6C  Ais*CT=(ABS(YMNl)-YN)/{XNNl-XNN) 

0  I  F.w  =  Atf  B-AMOT 
C I FY=YN-Y0 
DIFX=CIFY/DIFM 
XP=C I FX+XNN 
OIFFY=AMB*DIFX 
YP  =  C IFFY+YO 

IF  (YPNl.LT.O. )  YP=-YP 
IF  (ABSIXP-XND.LE.O.OQOl)  GO  TO  3501 
XN1 =XP 

WFSMljsWFSMBfl. 

IF  (WFSPfi.LE.25.  )  GO  TO  57 
GO  TO  3501 
35C0  CONTINUE 
XP^XHL 
YP---YFL 

3501  I F ( T { 5 )  -XNN)3412, 3412,3416 

3412  I F ( T ( 6 )  -  YN) 3413, 3413, 3416 

3 4 1 3  IF(ni)  -XNN) 34 16, 34  14, 3414 

3414  I F ( T ( 2 )  -  YNJ3416, 3415, 3415 

3415  JP  =  JP  ♦  1 

I P8  -  IPA  4-  JP 

XPL  r ( I PB )  =  XP 

YPLT( IPB)  =  YP 

IF  ( IPB. LT. 3099)  GO  TO  3416 

JP  =  0 

3  4  16  IC=ICH 

XPLCT ( IO)=XP 
YPLOT (10)  =  YP 

IF  (XNNl.LT.XTRAPl)  GO  TO  802 
IF  (XNNI.GT.XTRAP2JGC  TO  802 
IF  (ABS(YPNl).LT.YrRAPl)  GO  TO  802 
IF  (ARS(WN1  KGT.YTRAP2)  GO  TO  802 
WRITE  (6,660) 

GO  TC  102 

8C2  IF  (RESIL.GQ.O.O)  GO  TO  9781 
AL  PH  =  A  T  AN ( AMB  ) 

I3ETA=1, 57079 


I F { XNN1 . NE»XNN )  9ET4=ATAN (  AMDT) 

IF  { VCX'w  . GE . 0 .0  }  GO  TO  912 
IF  ( A3S { ALPHA-BE T4 )- 1 . 57079 }  810*904,308 
8C4  SIGNW.O 
S IGN2= 1 .0 
EPS=-1. 57079 
806  VB I T=0 . 

VOW  1  =  SORT { VDXW*VDXW\1+VDYW*VDYWN1) 

VB2N=VCV»  1  *RES  IL 
GO  rc  818 
808  S IG\i=l,0 
S IGN2  =  1,0 
GO  TO  816 
81C  SIGMs-l.o 
S I GN2--1 , 0 
GO  TG  816 

812  IF  U8S(4LPH-0ETA)-i. 57079)  808,814,810 
814  S IGN1= 1.0 
S I GN2= 1.0 
FPS=-l. 57079 
GO  TG  806 

816  VDW1  =  SGRT(A8S(  VCXW*  VCXWN 1 )  fABS  (  VDYw*  VOYWN 1  )  ) 
VRIN=V0W1*SIN( ALPH-9ET4) 

VBIT=VDW1*C0S(ALP! — BETA) 

V82\=V8IN»RESIL 
EPS  =  ATAK'(  VB2N/V8I  T  ) 

818  ANGLE=ALPH+EPS 

VDN2=SCRT{V8IT*VBIT+V8?N«V82N) 

VOX'XN  i=VCW2*S  IG\  1  *CQ$  (  ANGLE  ) 

V0YW.N1  =V0W2»S  ICN2*S  INI  ANGLE  ) 

IF  (YP.LT. 0.0)  V0YWN1=-VDYWN1 

PCu\CF= BOUNCE *1.0 

yl:ri:p=yp 

XL  =  XP 

V  YW  =  VDYWN I #u I nf 
VXi-<  =  VDXWNl»UlNF 

IF  (XL.GE.XMAX)  GO  TO  102 
IF  (BOUNCE. GT.  15.0)  GO  TO  102 
GO  TO  315 
9781  X  6  6  =  X  6  6 

IF  (XPF1.NE.0.0)  GO  TO  468 
IF  ( rN.EQ.O.O)  GG  TO  4  6  3 
FN=FNf 1,0 
GC  TO  464 

463  SAVE=YCROP 

r\=wo 

464  Y0PCP=YDR0P~CELYCP/2.0»»FN 
IF  (FN.LT.6.Q)  GO  TO  306 

Y  T  A\  =  YDROP 
YORCP=SAVE-CELYOP 
GO  f 0  1644 

468  IF  (FN.EQ.O.O)  GO  TO  61 
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FN=FN+  i  .0 

V 0 R C P = V f; R C P  <• C L  Y  C  P  /  2  -  0  «  *  F N 
IF  l  FI'S,  U  .6  »  }  GO  TO  306 

Y  T  AN= YOROP 

Y  C  P.  0  P  =  S  A  Y  £  -  0  £  L  Y  0  P 
GC  TO  61 

62  Y PV l = YP 
XPP1=XP 
UNC -UNDER 
WRITE  (6,3107) 

WRITE  (6, 124C)WATRAT,AV$,XP,YP, WATER 
WRITE  (6,3107) 

NK  =  NK  *■  1 

IF  (NK.EQ.U  GQ  TO  4  75 
l DX=4+2* ( NK- 1 ) 

TA310I IDX)=AVS 
T  AB 10  I IDX  +  I ) =WATRAT 
T  AB 10 ( 3 ) =NK 
XP=0.0 

475  IF  (FN. 60.0.0)  GC  TO  102 
YP=0 . 0 
XPP1=0.0 
GO  TC  6175 
61  Ah  I  T E  =  XP-XPM  1 
W  IC  =  YP-YPMl 
SGF I TE  =  AH  I TE«»2 
SOW  I C  =  W  1 0* *2 

AHYPT=SQRT( SGHIT6+SQW1D) 

IF(ALAMOA)  1641,  164C,  1641 
1640  CCN AER=AHYPT 

Z  =  60.0*  CELYOP 
IF  (FN.EG.O.O)  GC  TO  470 
Z=60,0»( YTAN-YCRLP) 

GO  TO  1642 

4/C  IF  (PPM. EG. 0.0)  GO  TO  1642 
Z =60.0* (YDROP-Y TAN) 

PPP  =  0 

GC  TO  1642 

164  1  SU"1AC  =  ABS(  {YP)MYPMU) 
CONAER=SUPRAC*AHYPT 
Z=AOc«  (ABS(2, »YCRCP-CELYDP  )  )*OELYOP 
IF  (FN.EO.O.O)  GO  TO  6165 
Z=60,0*{ YTAN-YCROP ) * ( YTAN+YDROP ) 

GO  TO  6167 

‘.-165  IF  ( PPM.  EQ.C.O)  GO  TO  6167 

Z  =  60. 0  * ( YDROP-Y  TAN ) * ( YOROP  hYT  AN ) 

PPM  =  0 . 0 

416/  IF  ( YP *YPM 1 . G£ . 0 „ 0 )  GQ  TO  1642 
PHYP  =  SQRT(YP*YP+f XP-XHL 1 )»*2  ) 
P"HYP=SQRT  (  YPM1*»2M  XPMI-XHH  )**2) 
CONAE  R= ( YP*PHYP ) ♦ ABS ( YPM 1 ■ PMHYP ) 

1642  WATER=Z« 12.0*UINF*ALWC 


WAFRAT=WATER/CONA£R 
164 4  XSS  =  XP-XHL 

XXSS  =  X PM 1  -  XHL 

CALL  TLUXXS  ( TAB3»XP*0»0,YBi) 

IF  ( ABS { YP J .LT. { Y31 f .001 } )  GO  TO  67 
CALL  TLUXXS  (TA86,XP,0,0fY32) 
YB2XP=YB2 

IF  { ABS(YP) .LT.(Y82XP+.001) }  GO  TO  68 
CALL  TLUXXS  { TA82 , X$ S , 0 , 0 , SUP ) 

UNOER  =- 1.0 

IF  (UNOERX.LT #0. }  GO  TO  64 
GO  TC  63 

68  CALL  TLUXXS  ( T A3 1 , XSS , 0 , 0 , SUP ) 
SUP=-$UP 
UND6R=1.0 

IF  (UNCERX.LT. O.C)  C-O  TO  64 

63  CALL  TLUXXS ( TAB l , XXS S , 0 , 0 »  SUPM 1 ) 

SUP*V  1  =-SUPM  1 

GO  TO  65 

64  CALL  TLUXXS  ( TAB2 , X X SS . 0 , 0 , SUPM 1 ) 

65  AVS  =  (SUP  +  SUPMIJ/2.0 
IF  ( XPM1 .NE.0.0 )  GO  TO  62 

6175  AVS  =  SUP 

PMv=1>0 

WATRAT=0,0 
FN  = 0.0 
GC  TO  62 

67  CALL  TLUXXS  ( T AS9 , X S S , 0 , 0, SUP ) 

CALL  TLUXXS  (  TA89,XXSSf 0»C, SUPMl ) 

IF  (YP.GT.O.)  GO  TO  70 
UNDER =-t. 

SUP=-SUP 
GC  TC  71 
fZ  U  N  0  £  R  =  I . 

?I  IF  (UNCERX.LT. 0.0)  SUPM 1 =- SUPM 1 
GO  TC  65 
END 
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